******************************************************************************************************************************************************************************
* Do file for
* Do Physically Stronger Males Prevail in Non-Physical Conflicts?
* By Dan Nguyen, Michael Bang Petersen, Julia Nafziger, and Alexander K. Koch
* Aarhus University
* June 2020
* Contact: Alexander K. Koch, akoch@econ.au.dk
******************************************************************************************************************************************************************************

**********************************************************
***********Study 1***********
**********************************************************
use data_study1.dta, clear

*note:
*study included rating of the following types pictures of 50 targets:
*- body silhouettes (BS)
*- full body (FB)
*- full face (FF)
*- face silhouette (FS)
*200 raters rated body pictures, both full and silhouettes in 2 sets with random ordering within each set; 100 raters started with FB and then did FS, 100 had reverse order of sets
*200 raters rated face pictures, both full and silhouettes in 2 sets with random ordering within each set; 100 raters started with FB and then did FS, 100 had reverse order of sets

reshape long StrFS@ StrFF@ StrBS@ StrFB@ AttFS@ AttFF@ AttBS@ AttFB@ age_@ bicep_@ fighting_@ height_@ maxchest_@ maxgrip_@ perceivedattractiveness_@ perceivedstrength_@ semester_@ strengthtrain_@ study_@ weight_@, i(Participant) j(rated)


*standardize variables: 
*fix row in data and take mean and std.dev. over all the rated persons-> gives same result for each row 
by Participant, sort: egen M=mean(bicep_)
by Participant, sort: egen SD=sd(bicep_)
by rated, sort: gen zbicep= (bicep_ - M)/SD
drop M SD

by Participant, sort: egen M=mean(height_)
by Participant, sort: egen SD=sd(height_)
by rated, sort: gen zheight = (height_ - M)/SD
drop M SD

by Participant, sort: egen M=mean(maxchest_)
by Participant, sort: egen SD=sd(maxchest_)
by rated, sort: gen zmaxchest = (maxchest_ - M)/SD
drop M SD

by Participant, sort: egen M=mean(maxgrip_)
by Participant, sort: egen SD=sd(maxgrip_)
by rated, sort: gen zmaxgrip = (maxgrip_ - M)/SD
drop M SD


by Participant, sort: egen M=mean(perceivedattractiveness_)
by Participant, sort: egen SD=sd(perceivedattractiveness_)
by rated, sort: gen zperceivedatt = (perceivedattractiveness_ - M)/SD
drop M SD

by Participant, sort: egen M=mean(perceivedstrength_)
by Participant, sort: egen SD=sd(perceivedstrength_)
by rated, sort: gen zperceivedstr = (perceivedstrength_ - M)/SD
drop M SD

** standardized measure for composite strength based on Sell et al.**
gen zform = ( zbicep+ zmaxchest+ zmaxgrip+ zperceivedstr)/4

** We are only interested in body silhouette pictures and strength ratings, drop 200 raters who only rated faces**
drop if StrBS==.

*summary stats on age
*raters of body silhouettes (n=200)
by Participant, sort: gen particobs=_n
sum age if particobs==1
*rated targets (n=50)
by rated, sort: gen ratedobs=_n
sum age_ if ratedobs==1
*descriptive stats for measurements of targets (ESM)
sum height_ weight_ bicep_ maxchest_ perceivedstrength_ if ratedobs==1

** Cronbach Alpha for composite strength (ESM) **
alpha zbicep zmaxchest zmaxgrip zperceivedstr, item

*test physical strength in group A vs B
ttest zform if ratedobs==1, by(Treatment) 


**Correlations** Replication of Sell
by rated, sort: egen BSmean=mean(StrBS)
by Participant, sort: egen M=mean(BSmean)
by Participant, sort: egen SD=sd(BSmean)
*z-score mean ratings per rated
by rated, sort: gen zBSmean = (BSmean - M)/SD
drop M SD
*z-score individual ratings 
egen M=mean(StrBS)
egen SD=sd(StrBS)
gen zStrBS = (StrBS - M)/SD
drop M SD
*correlation of average rating of strength from body silhouettes with the actual physical strength measure
pwcorr zform BSmean if ratedobs==1, sig ob // include each of the 50 rated individuals only once, i.e. 50 observations
*average individual accuracy: coefficient from regression of actual strength measure on rating, cluster at rater (200 raters) and rated (50 rated) level 
vce2way reg zform zStrBS, cluster(Participant rated)


******************************************************************************************************************************************************************************
******************************************************************************************************************************************************************************

**********************************************************
***********Study 2***********
**********************************************************

use data_study2,clear
set more off

cd "./results"
**********************************************************
***********Prepare estimation ***********
**********************************************************

which vce2way //check if user written program is installed, if not use findit <program name> in stata 
*does not allow if option, hence below we need to temporarily drop all observations not in estimation sample 

*For tables generate labels for interaction variables
label var absdstrength "|Difference in strength|"
label var absdheight "|Difference in height|"
label var absdstrength2 "|Difference in strength|^2"

*figure scheme
set scheme lean1

********************************************************************************
*********** Stats/results in manuscript text ***********
********************************************************************************

sum age if period==1

** Cronbach Alpha for composite strength **
alpha z_flexedbicep z_maxgrip z_maxchest z_perceivedstrength if period==1, item

** Conflict duration **
sum finalcosts if treatment==1 & contestantincluded==1 //MAT
sum finalcosts if treatment==0 & contestantincluded==1 //SAT

** Regression coefficients reported / marginal effects **
*DURATION
quiet: sum finalcost if treatment==1 & contestantincluded==1
local sdhat1=r(sd)
quiet: sum absdstrength if contestantincluded==1
local sdhat2=r(sd)
preserve 
keep if treatment==1 & contestantincluded==1 
vce2way nbreg finalcosts absdstrength  per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 , cluster(id1 id2)
** standardized average marginal effect **
margins, dydx(*) post
matrix temp=e(b)
matrix list temp
local bhat=temp[1,1]
local effectsize=(`bhat'*`sdhat2')/`sdhat1'
display `effectsize' "(standardized average marginal effect)"
restore

*WINNER
*Test that absdstrength and absdstrength2 are jointly 0 in regression for winner*
preserve
keep if treatment==1 & contestantincluded==1 & draw!=1
quietly: vce2way logit strongwins absdstrength c.absdstrength#c.absdstrength per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 , cluster(id1 id2)
test absdstrength c.absdstrength#c.absdstrength
restore
*Discrete marginal effect on prob(strongest wins)
*percentage point increase in the probability of the stronger contestant winning when the difference in strength increases by one standard deviation from its mean
preserve
quiet: sum absdstrength if contestantincluded==1
local sdhat2=r(sd)
keep if treatment==1 & contestantincluded==1 
vce2way logit strongwins absdstrength c.absdstrength#c.absdstrength  per1 per2 per3 per4 per5 per6 per7 per8 per9 per10, cluster(id1 id2)
mtable, atmeans at(absdstrength = gen(absdstrength))  ///
		at(absdstrength = gen(absdstrength+`sdhat2'))  post
mlincom 2 - 1, rowname(1std.dev.increase |diff.strength|)
matrix temp=_mlincom
matrix list temp
local effectsize=temp[1,1]
display `effectsize' "(average change in probability that stronger contestant wins for one std.dev. change in |difference strength|)"
restore
*Discrete marginal effect on prob(strongest wins), zero duration contests
preserve
quiet: sum absdstrength if contestantincluded==1
local mean2=r(mean)
local sdhat2=r(sd)
local mean2plus= `mean2'+`sdhat2'
keep if treatment==1 & contestantincluded==1 & finalcost==0
vce2way logit strongwins absdstrength per1 per2 per3 per4 per5 per6 per7 per8 per9 per10, cluster(id1 id2)
mtable, atmeans at(absdstrength = gen(absdstrength))  ///
		at(absdstrength = gen(absdstrength+`sdhat2'))  post
mlincom 2 - 1, rowname(1std.dev.increase |diff.strength|)
matrix temp=_mlincom
matrix list temp
local effectsize=temp[1,1]
display `effectsize' "(average change in probability that stronger contestant wins for one std.dev. change in |difference strength|; zero duration contests)"
restore





********************************************************************************
***********Fig. 2 (left panel): The effect of differences in strength on conflict duration in the Mutual Assessment Treatment***********
********************************************************************************
*lines in plot area for 5th and 95th percentile of absdstrength
sum absdstrength
centile absdstrength, centile(5 95)

preserve 
keep if treatment==1 & contestantincluded==1 
quietly: vce2way nbreg finalcosts absdstrength  per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 , cluster(id1 id2)
quietly: margins, at(absdstrength=(0(0.1)3.9)) atmeans
marginsplot,xline(.0539063) xline( 1.910426) recast(line) recastci(rline) plotopt(lw(medthick) lcolor(dknavy)) ciopt(lpat(dash) lw(vthin)) ytitle("Predicted duration in seconds") xtitle("Absolute difference in strength") title("", size(small)) xscale(titlegap(3)) ylabel(0(20)80) yscale(r(-5,90) titlegap(2))  saving(durationSTRonly, replace) 
restore


********************************************************************************
***********Fig. 3 (left panels): Who wins the contests in the Mutual Assessment Treatment***********
********************************************************************************
*All contests
preserve
keep if treatment==1 & contestantincluded==1 & draw!=1
quietly: vce2way logit strongwins absdstrength c.absdstrength#c.absdstrength per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 , cluster(id1 id2)
quietly: margins, at(absdstrength=(0(0.1)3.9)) atmeans
marginsplot,xline(.0539063) xline( 1.910426) recast(line) recastci(rline) plotopt(lw(medthick) lcolor(dknavy)) ciopt(lpat(dash) lw(vthin)) ytitle("Probability of stronger winning") xtitle("Absolute difference in strength") title("a) Strongest wins", size(small))  xscale(titlegap(3)) ylabel(0(.2)1) yscale(r(-.3,1) titlegap(2)) scheme(lean1) saving(winner,replace)  
restore
*zero duration contests
preserve
keep if treatment==1 & contestantincluded==1 & finalcost==0
quietly: vce2way logit strongwins absdstrength per1 per2 per3 per4 per5 per6 per7 per8 per9 per10, cluster(id1 id2)
quietly: margins, at(absdstrength=(0(0.1)3.9)) atmeans
marginsplot,xline(.0539063) xline( 1.910426) recast(line) recastci(rline) plotopt(lw(medthick) lcolor(dknavy)) ciopt(lpat(dash) lw(vthin)) ytitle("Probability of stronger winning") xtitle("Absolute difference in strength") title("c) Strongest wins (zero duration contests)", size(small))  xscale(titlegap(3)) ylabel(0(.2)1) yscale(r(-.3,1.2) titlegap(2)) scheme(lean1) saving(winnerzero,replace)  
restore


***********************************************
***********Fig. 4 (except right panels): Evidence for mutual assessment***********
***********************************************
gen dummyregressor1=finalcosts
gen dummyregressor2=-finalcosts

*Joint test of both coefficients being zero
preserve 
keep if treatment==1 & contestantincluded==1 
vce2way nbreg finalcosts strongformid weakformid per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 , cluster(id1 id2)
test strongformid_z  weakformid_z
restore

preserve 
keep if treatment==1 & contestantincluded==1 
quiet: vce2way nbreg finalcosts strongformid weakformid per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 , cluster(id1 id2)
quiet: margins, at(weakformid=(-1.61(0.1)0.831)) atmeans
quiet: marginsplot, recast(line) recastci(rline) ciopts(lpattern(dash)) plotopt(lcolor(dknavy) lw(thick)) xlabel(#5) legend(off) scheme(lean1) ytitle("Predicted duration in seconds") xtitle("Strength of weaker contestant" "(standardized)") title("c) Observed results in Study 2" " ") xscale(titlegap(3)) yscale(titlegap(2)) saving(weakformid,replace)
quiet: vce2way nbreg finalcosts strongformid weakformid per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 , cluster(id1 id2)
quiet: margins, at(strongformid=(-1.26(0.1)2.30)) atmeans
quietly: marginsplot, recast(line) recastci(rline) ciopts(lpattern(dash)) plotopt(lcolor(dknavy) lw(thick)) xlabel(#5) legend(off) scheme(lean1) ytitle("Predicted duration in seconds") xtitle("Strength of stronger contestant" "(standardized)") title("") saving(strongformid,replace)
graph twoway (function x/2, range(strongform) n(2) lpattern(dash)  lpattern(dash) lcolor(dknavy) lw(vthick) xlabel(#5) legend(off) scheme(lean1) ytitle("Duration") yscale(range(-3 3))  ylabel(none ) xlabel(none -1.2"Low" 2.2"High") title("a) Theoretical predictions:" "Mutual Assessment Model") xtitle("Strength of weaker contestant" "") saving(weak_dummy1,replace))
graph twoway (function -x/2, range(strongform) n(2) lpattern(dash)  lpattern(dash) lcolor(dknavy) lw(vthick) xlabel(#5) legend(off) scheme(lean1) ytitle("Duration") yscale(range(-3 3))  ylabel(none ) xlabel(none -1.2"Low" 2.2"High") title("") xtitle("Strength of stronger contestant" "") saving(strong_dummy1,replace))
graph twoway (function x/2, range(strongform) n(2) lpattern(dash)  lpattern(dash) lcolor(dknavy) lw(vthick) xlabel(#5) legend(off) scheme(lean1) ytitle("Duration") yscale(range(-3 3))  ylabel(none ) xlabel(none -1.2"Low" 2.2"High") title("b) Theoretical predictions:" "Pure Self-Assessment Model") xtitle("Strength of weaker contestant" "") saving(weak_dummy2,replace))
graph twoway (function x/4, range(strongform) n(2) lpattern(dash)  lpattern(dash) lcolor(dknavy) lw(vthick) xlabel(#5) legend(off) scheme(lean1) ytitle("Duration") yscale(range(-3 3))  ylabel(none ) xlabel(none -1.2"Low" 2.2"High") title("") xtitle("Strength of stronger contestant" "") saving(strong_dummy2,replace))
restore




***************************************************************
***********Figure S3 Study 2: Distribution of the marginal effects ***********
***************************************************************
*distribution of marginal effects
sum finalcost if treatment==1 & contestantincluded==1
local sdhat1=r(sd)
sum absdstrength if contestantincluded==1
local sdhat2=r(sd)
preserve 
keep if treatment==1 & contestantincluded==1 
vce2way nbreg finalcosts absdstrength  per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 , cluster(id1 id2)
*standardized marginal effect 
margins, generate(MEduration) ///
at(absdstrength = gen(absdstrength)) at(absdstrength = gen(absdstrength+`sdhat2'))
gen DCduration = (MEduration2 - MEduration1)/`sdhat1'
hist DCduration, ytitle("Density") xtitle("Standardized marginal effect on contest duration") title("a) Contest duration") xscale(titlegap(3)) yscale(titlegap(2)) scheme(lean1) saving(MEduration, replace) 
sum DCduration
restore

*distribution of increase in prob(stronger wins) for one std.dev. increase in |diff strength|
preserve
quiet: sum absdstrength if contestantincluded==1
local sdhat2=r(sd)
keep if treatment==1 & contestantincluded==1 & draw!=1
vce2way logit strongwins absdstrength  c.absdstrength#c.absdstrength per1 per2 per3 per4 per5 per6 per7 per8 per9 per10, cluster(id1 id2)
*standardized marginal effect at means of covariates
margins, generate(MEwinner) ///
at(absdstrength = gen(absdstrength)) at(absdstrength = gen(absdstrength+`sdhat2'))
gen DCwinner = MEwinner2 - MEwinner1
hist DCwinner, ytitle("Density") xtitle("Change in prob(strongest wins)") title("b) Strongest wins") xscale(titlegap(3)) yscale(titlegap(2)) scheme(lean1) saving(MEwinner, replace) 
restore

*distribution of increase in prob(stronger wins) for one std.dev. increase in |diff strength|; for zero duration contests
preserve
quiet: sum absdstrength if contestantincluded==1
local sdhat2=r(sd)
keep if treatment==1 & contestantincluded==1 & finalcost==0
vce2way logit strongwins absdstrength per1 per2 per3 per4 per5 per6 per7 per8 per9 per10, cluster(id1 id2)
*standardized marginal effect at means of covariates
margins, generate(MEwinnerzero) ///
at(absdstrength = gen(absdstrength)) at(absdstrength = gen(absdstrength+`sdhat2'))
gen DCwinnerzero = MEwinnerzero2 - MEwinnerzero1
hist DCwinnerzero, ytitle("Density") xtitle("Change in prob(strongest wins)") title("c) Strongest wins (zero duration contests)") xscale(titlegap(3)) yscale(titlegap(2)) scheme(lean1) saving(MEwinnerzero, replace) 
restore

*Distribution of the effects of a one standard deviation" "increase in the absolute difference in strength
gr combine MEduration.gph MEwinner.gph MEwinnerzero.gph, ///
	col(2) hole(2) iscale(0.7) xcommon title("") scheme(lean1) 
graph export figureS3.pdf, fontface(Arial) replace





**********************************************************
***********Results (ESM) ***********
**********************************************************

*********** Number of contests (ESM) ***********
* tabulate the maximum number of periods against sessions: number of participants=number of periods+1
by session period, sort: gen sessnum=_N
tab sessnum if period==1 & treatment==1  //MAT
tab sessnum if period==1 & treatment==1  //SAT

*********** Coefficients on the stronger contestant and weaker contestants’ strength (ESM) ***********
*the coefficients on the stronger contestant and weaker contestants’ strength have the same absolute value*

preserve 
*model 7
keep if treatment==1 & contestantincluded==1
vce2way nbreg finalcosts strongformid weakformid per1 per2 per3 per4 per5 per6 per7 per8 per9 per10   , cluster(id1 id2)
test strongformid+weakformid=0
restore


***********Correlation of Formidability and risk tolerance (ESM) ***********

*"we find no correlation between strength and risk tolerance"
preserve
collapse age risk maxgrip maxchest flexedbicep perceivedstrength perceivedattractiveness height weight Formidability_z bmi treatment session, by(participantid)
pwcorr  risk Formidability_z, sig obs
restore 

**********************************************************
***********Table S1: Descriptive statistics for Study 2 (non-standardized) ***********
**********************************************************

*Mutual Assessment Treatment
sum  Formidability_z maxgrip maxchest flexedbicep perceivedstrength age height weight bmi risk if treatment==1 & period==1
*Self-Assessment Treatment
sum  Formidability_z maxgrip maxchest flexedbicep perceivedstrength age height weight bmi risk if treatment==0 & period==1

*Balance tests and CIs
preserve
collapse age risk maxgrip maxchest flexedbicep perceivedstrength perceivedattractiveness height weight Formidability_z bmi treatment session, by(participantid)
ttest age , by(treatment)
ttest Formidability_z , by(treatment)
ttest height , by(treatment)
ttest weight , by(treatment)
ttest bmi , by(treatment)
ttest maxgrip , by(treatment)
ttest maxchest , by(treatment)
ttest flexedbicep , by(treatment)
ttest perceivedstrength , by(treatment)
ttest risk , by(treatment)
restore 

**********************************************************
***********Table S2: Negative binomial regressions for the duration of contests (Mutual Assessment Treatment in Study 2) ***********
**********************************************************
which regsave //user written stata command for saving regression output
*Prepare a matrix for 99-percent CI*
set more off
matrix alphaCI = J(4,8,.)
matrix colnames alphaCI =  model1 model2 model3 model4 model5 model6 model7 model8
matrix rownames alphaCI = alpha alphase lower upper

preserve 
keep if treatment==1 & contestantincluded==1
set level 99
*model1
vce2way nbreg finalcosts absdstrength per1 per2 per3 per4 per5 per6 per7 per8 per9 per10   , cluster(id1 id2)
eststo model1
matrix tempmat=r(table)
local colno=colsof(tempmat)
matrix alphaCI[1,1]=tempmat[1,`colno']
matrix alphaCI[2,1]=tempmat[2,`colno']
matrix alphaCI[3,1]=tempmat[5,`colno']
matrix alphaCI[4,1]=tempmat[6,`colno']
*model2
vce2way nbreg finalcosts absdheight per1 per2 per3 per4 per5 per6 per7 per8 per9 per10   , cluster(id1 id2)
eststo model2
matrix tempmat=r(table)
local colno=colsof(tempmat)
matrix alphaCI[1,2]=tempmat[1,`colno']
matrix alphaCI[2,2]=tempmat[2,`colno']
matrix alphaCI[3,2]=tempmat[5,`colno']
matrix alphaCI[4,2]=tempmat[6,`colno']
*model3
vce2way nbreg finalcosts absdstrength absdheight c.absdstrength#c.absdheight per1 per2 per3 per4 per5 per6 per7 per8 per9 per10   , cluster(id1 id2)
eststo model3
matrix tempmat=r(table)
local colno=colsof(tempmat)
matrix alphaCI[1,3]=tempmat[1,`colno']
matrix alphaCI[2,3]=tempmat[2,`colno']
matrix alphaCI[3,3]=tempmat[5,`colno']
matrix alphaCI[4,3]=tempmat[6,`colno']
*model4
vce2way nbreg finalcosts absdbmi per1 per2 per3 per4 per5 per6 per7 per8 per9 per10   , cluster(id1 id2)
eststo model4
matrix tempmat=r(table)
local colno=colsof(tempmat)
matrix alphaCI[1,4]=tempmat[1,`colno']
matrix alphaCI[2,4]=tempmat[2,`colno']
matrix alphaCI[3,4]=tempmat[5,`colno']
matrix alphaCI[4,4]=tempmat[6,`colno']
*model5
vce2way nbreg finalcosts absdrisk per1 per2 per3 per4 per5 per6 per7 per8 per9 per10   , cluster(id1 id2)
eststo model5
matrix tempmat=r(table)
local colno=colsof(tempmat)
matrix alphaCI[1,5]=tempmat[1,`colno']
matrix alphaCI[2,5]=tempmat[2,`colno']
matrix alphaCI[3,5]=tempmat[5,`colno']
matrix alphaCI[4,5]=tempmat[6,`colno']
*model6
vce2way nbreg finalcosts absdstrength absdrisk per1 per2 per3 per4 per5 per6 per7 per8 per9 per10   , cluster(id1 id2)
eststo model6
matrix tempmat=r(table)
local colno=colsof(tempmat)
matrix alphaCI[1,6]=tempmat[1,`colno']
matrix alphaCI[2,6]=tempmat[2,`colno']
matrix alphaCI[3,6]=tempmat[5,`colno']
matrix alphaCI[4,6]=tempmat[6,`colno']
**model7: leave out as minor detail beyond models 5 and 8
*vce2way nbreg finalcosts strongrisk weakrisk per1 per2 per3 per4 per5 per6 per7 per8 per9 per10   , cluster(id1 id2)
*eststo model7
*matrix tempmat=r(table)
*local colno=colsof(tempmat)
*matrix alphaCI[1,7]=tempmat[1,`colno']
*matrix alphaCI[2,7]=tempmat[2,`colno']
*matrix alphaCI[3,7]=tempmat[5,`colno']
*matrix alphaCI[4,7]=tempmat[6,`colno']
*model7
vce2way nbreg finalcosts strongformid weakformid per1 per2 per3 per4 per5 per6 per7 per8 per9 per10   , cluster(id1 id2)
eststo model7
matrix tempmat=r(table)
local colno=colsof(tempmat)
matrix alphaCI[1,7]=tempmat[1,`colno']
matrix alphaCI[2,7]=tempmat[2,`colno']
matrix alphaCI[3,7]=tempmat[5,`colno']
matrix alphaCI[4,7]=tempmat[6,`colno']
*model8
vce2way nbreg finalcosts strongrisk weakrisk strongformid weakformid taller_z shorter_z per1 per2 per3 per4 per5 per6 per7 per8 per9 per10   , cluster(id1 id2)
eststo model8
matrix tempmat=r(table)
local colno=colsof(tempmat)
matrix alphaCI[1,8]=tempmat[1,`colno']
matrix alphaCI[2,8]=tempmat[2,`colno']
matrix alphaCI[3,8]=tempmat[5,`colno']
matrix alphaCI[4,8]=tempmat[6,`colno']
restore
matrix list alphaCI
esttab using duration_MAT.rtf, replace /// 
	addnote("Robust standard errors") drop(per*) label nonumber mtitle("Model 1" "Model 2" "Model 3" "Model 4" "Model 5" "Model 6" "Model 7" "Model 8") b(3) se(3) se starlevels(* 0.1 ** 0.05 *** 0.01)  
*note: 38 clusters because there are 4 sessions and the max and min subject ids from each session only appear once, either in id1 or id2, whereas all other subject ids appear in both id1 and id2
eststo clear

**********************************************************
*model selection
**********************************************************
*Note: Can't use cluster with LR test as the likelihood does not account for the clustering, alternative is to use Wald test, 
*ie. CI reported for alpha, level set to 1% (see above)

*LR test
*model1
poisson finalcosts absdstrength per1 per2 per3 per4 per5 per6 per7 per8 per9 per10  if treatment==1 & contestantincluded==1 , robust
est store poisson
nbreg finalcosts absdstrength per1 per2 per3 per4 per5 per6 per7 per8 per9 per10  if treatment==1 & contestantincluded==1 , robust
est store nbreg
lrtest poisson nbreg, stats force
eststo clear

*model2
poisson finalcosts absdheight per1 per2 per3 per4 per5 per6 per7 per8 per9 per10   if treatment==1 & contestantincluded==1, robust
est store poisson
nbreg finalcosts absdheight per1 per2 per3 per4 per5 per6 per7 per8 per9 per10   if treatment==1 & contestantincluded==1, robust
est store nbreg
lrtest poisson nbreg, stats force
eststo clear

*model3
poisson finalcosts absdstrength absdheight c.absdstrength#c.absdheight per1 per2 per3 per4 per5 per6 per7 per8 per9 per10   if treatment==1 & contestantincluded==1, robust
est store poisson
nbreg finalcosts absdstrength absdheight c.absdstrength#c.absdheight per1 per2 per3 per4 per5 per6 per7 per8 per9 per10   if treatment==1 & contestantincluded==1, robust
est store nbreg
lrtest poisson nbreg, stats force
eststo clear

*model4
poisson finalcosts absdbmi per1 per2 per3 per4 per5 per6 per7 per8 per9 per10  if treatment==1 & contestantincluded==1 , robust
est store poisson
nbreg finalcosts absdbmi per1 per2 per3 per4 per5 per6 per7 per8 per9 per10  if treatment==1 & contestantincluded==1 , robust
est store nbreg
lrtest poisson nbreg, stats force
eststo clear


*model5
poisson finalcosts absdrisk per1 per2 per3 per4 per5 per6 per7 per8 per9 per10   if treatment==1 & contestantincluded==1, robust
est store poisson
nbreg finalcosts absdrisk per1 per2 per3 per4 per5 per6 per7 per8 per9 per10   if treatment==1 & contestantincluded==1, robust
est store nbreg
lrtest poisson nbreg, stats force
eststo clear

*model6
poisson finalcosts absdstrength absdrisk per1 per2 per3 per4 per5 per6 per7 per8 per9 per10   if treatment==1 & contestantincluded==1, robust
est store poisson
nbreg finalcosts absdstrength absdrisk per1 per2 per3 per4 per5 per6 per7 per8 per9 per10   if treatment==1 & contestantincluded==1, robust
est store nbreg
lrtest poisson nbreg, stats force
eststo clear

*model7
poisson finalcosts strongformid weakformid per1 per2 per3 per4 per5 per6 per7 per8 per9 per10   if treatment==1 & contestantincluded==1, robust
est store poisson
nbreg finalcosts strongformid weakformid per1 per2 per3 per4 per5 per6 per7 per8 per9 per10   if treatment==1 & contestantincluded==1, robust
est store nbreg
lrtest poisson nbreg, stats force
eststo clear

*model8
poisson finalcosts strongrisk weakrisk strongformid weakformid taller_z shorter_z per1 per2 per3 per4 per5 per6 per7 per8 per9 per10   if treatment==1 & contestantincluded==1, robust
est store poisson
nbreg finalcosts strongrisk weakrisk strongformid weakformid taller_z shorter_z per1 per2 per3 per4 per5 per6 per7 per8 per9 per10   if treatment==1 & contestantincluded==1, robust
est store nbreg
lrtest poisson nbreg, stats force
eststo clear

*CONCLUSION: alpha significantly differs from 0 -> overdispersion, reject poisson in favor of neg.binomial model


************************************************************************
***********Table S4: Logistic regressions on whether the strongest contestant wins (Mutual Assessment Treatment in Study 2)**********
************************************************************************


preserve
keep if treatment==1 & contestantincluded==1 & draw!=1 
vce2way logit strongwins absdstrength per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 ,  cluster(id1 id2)
eststo model1
vce2way logit strongwins absdstrength absdstrength2 per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 , cluster(id1 id2)
eststo model2
//To test that absdstrength and absdstrength are jointly 0:
test absdstrength absdstrength2

vce2way logit strongwins strongrisk weakrisk per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 , cluster(id1 id2)
eststo model3
vce2way logit strongwins absdstrength absdstrength2 strongrisk weakrisk per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 , cluster(id1 id2)
eststo model4
keep if treatment==1 & contestantincluded==1 & finalcosts==0 
vce2way logit strongwins absdstrength per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 ,  cluster(id1 id2)
eststo model5
vce2way logit strongwins strongrisk weakrisk per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 , cluster(id1 id2)
eststo model6
vce2way logit strongwins absdstrength strongrisk weakrisk per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 , cluster(id1 id2)
eststo model7
restore

esttab using winner_strength_MAT.rtf, replace ///
	addnote("Logistic regression") drop(per*) label nonumber mtitle("Model 1" "Model 2" "Model 3" "Model 4") b(3) se(3) se starlevels(* 0.10 ** 0.05 *** 0.01)
eststo clear


**********************************************************
***********Table S6: Contest outcomes Self-Assessment Treatment vs Mutual Assessment Treatment in Study 2.***********
**********************************************************
gen MAT=(treatment==1)

*as is standard practice, add main effects in addition to interaction effect

preserve
keep if (treatment==1 |treatment==0) & contestantincluded==1

vce2way nbreg finalcosts MAT i.MAT#c.absdstrength per1 per2 per3 per4 per5 per6 per7 per8 per9 per10  , cluster(id1 id2)
eststo model1
vce2way nbreg finalcosts MAT i.MAT#c.strongformid i.MAT#c.weakformid per1 per2 per3 per4 per5 per6 per7 per8 per9 per10  , cluster(id1 id2)
eststo model2
vce2way logit strongwins MAT i.MAT#c.absdstrength per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 , cluster(id1 id2)
eststo model3
vce2way logit strongwins MAT i.MAT#c.absdstrength i.MAT#c.absdstrength2 per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 , cluster(id1 id2)
eststo model4
restore 
esttab using MATSAT.rtf, pr2 replace /// 
	addnote("Robust standard errors") drop(per*) label nonumber mtitle("Model 1" "Model 2" "Model 3" "Model 4" "Model 5" "Model 6" "Model 7") b(3) se(3) se starlevels(* 0.1 ** 0.05 *** 0.01)  

eststo clear



***********************************************************************
***********Table S7: (Study 2) Logistic regressions on whether the tallest contestant wins (Mutual Assessment Treatment)****************
***********************************************************************

*exclude cases where sameheight==1
by gameid, sort: egen sumtallest = sum(tallest) 

preserve
keep if tallest==1 & sumtallest==1 & treatment==1 & draw!=1

vce2way logit tallwin absdheight per1 per2 per3 per4 per5 per6 per7 per8 per9 per10  , cluster(id1 id2)
eststo model1
vce2way logit tallwin absdheight absdheight2 per1 per2 per3 per4 per5 per6 per7 per8 per9 per10  , cluster(id1 id2)
eststo model2
vce2way logit tallwin tallrisk shortrisk per1 per2 per3 per4 per5 per6 per7 per8 per9 per10  , cluster(id1 id2)
eststo model3
vce2way logit tallwin absdheight tallrisk shortrisk per1 per2 per3 per4 per5 per6 per7 per8 per9 per10  , cluster(id1 id2)
eststo model4
restore 
esttab using winner_height_MAT.rtf, replace ///
	addnote("Logistic regression") drop(per*) label nonumber mtitle("Model 1" "Model 2" "Model 3" "Model 4") b(3) se(3) se starlevels(* 0.10 ** 0.05 *** 0.01) 
eststo clear

************************************************************************
***********Table S8: (Study 2) Logistic regressions on whether the leanest contestant wins (Mutual Assessment Treatment)********************
************************************************************************

*exclude cases where samebmi==1
by gameid, sort: egen sumleanest = sum(leanest)

preserve
keep if treatment==1 & leanest==1 & sumleanest==1 & draw!=1

vce2way logit leanwin absdbmi per1 per2 per3 per4 per5 per6 per7 per8 per9 per10  , cluster(id1 id2)
eststo model1
vce2way logit leanwin absdbmi absdbmi2 per1 per2 per3 per4 per5 per6 per7 per8 per9 per10  , cluster(id1 id2)
eststo model2
vce2way logit leanwin leanrisk heavyrisk per1 per2 per3 per4 per5 per6 per7 per8 per9 per10  , cluster(id1 id2)
eststo model3
vce2way logit leanwin absdbmi leanrisk heavyrisk per1 per2 per3 per4 per5 per6 per7 per8 per9 per10  , cluster(id1 id2)
eststo model4
restore 
esttab using winner_bmi_MAT.rtf,  replace ///
	addnote("Logistic regression") drop(per*) label nonumber mtitle("Model 1" "Model 2" "Model 3" "Model 4") b(3) se(3) se starlevels(* 0.10 ** 0.05 *** 0.01) 
eststo clear


******************************************************************************************************************************************************************************
******************************************************************************************************************************************************************************
**********************************************************
***********Power analysis for Study 3 (uses results from Study 2)***********
**********************************************************

*Note: we use pre-run Monte Carlo simulation


*Baseline from Study 2 (to get model coefficients; study 4 data set has all relevant variables defined)*

cd "..\"	
	 
use data_study4.dta, clear
set more off


keep if contestantincluded==1
vce2way nbreg finalcosts absdstrength  per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 , cluster(id1 id2)

*ICC
capture drop residact predact
				predict double predact, n
				gen residact=finalcosts-predact

loneway residact id1
scalar actirr=r(rho)
loneway residact id2
scalar actirr=(actirr+r(rho))/2
scalar list actirr
*R^2*
reg finalcosts absdstrength  per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 
reg finalcosts absdstrength  


**************************************************************************************
*Figure S5: Power as a function of coefficient size and sample size
**************************************************************************************
***With pre-run Monte Carlo results (generated using programs below)***
quiet{
quiet{
use mymontecarlolownoise, clear
matrix mcarloresults90 = J(11,3,.) 
matrix mcarloresults95 = J(11,3,.) 
matrix mcarloresults99 = J(11,3,.) 
 
local mcol=1
local coefvals ="03 04 05"
foreach x of local coefvals  {
	local mrow=1
	forvalues i=5(1)15 {
	 quiet sum sim`x'90_`i'
	 matrix mcarloresults90[`mrow',`mcol']=r(mean) 
	 quiet sum sim`x'95_`i'
	 matrix mcarloresults95[`mrow',`mcol']=r(mean) 
	 quiet sum sim`x'99_`i'
	 matrix mcarloresults99[`mrow',`mcol']=r(mean) 
	 local mrow=`mrow'+1
	 }
	local mcol=`mcol'+1
 }
 }
sum simest03
sum simr203
gen avicc03= (simicc103+simicc203)/2
sum avicc03

sum simest04
sum simr204
gen avicc04= (simicc104+simicc204)/2
sum avicc04

sum simest05
sum simr205
gen avicc05= (simicc105+simicc205)/2
sum avicc05
 
matrix nobs=(50, 60, 70, 80, 90, 100, 110, 120, 130, 140, 150)'
matrix ncontests=(225, 270, 315, 360, 405, 450, 495, 540, 585, 630, 675)'
matrix A==J(11,3,1)
matrix power90=A -mcarloresults90
matrix power95=A -mcarloresults95
matrix power99=A -mcarloresults99
matrix lownoise=(nobs,ncontests, power90, power95 ,power99)
matrix colnames lownoise =  N Contests .3p90 .4p90 .5p90 .3p95 .4p95 .5p95 .3p99 .4p99 .5p99
matrix list lownoise

** high noise **
quiet{
use mymontecarlohighnoise, clear
matrix mcarloresults90 = J(11,3,.) 
matrix mcarloresults95 = J(11,3,.) 
matrix mcarloresults99 = J(11,3,.) 
 
local mcol=1
local coefvals ="03 04 05"
foreach x of local coefvals  {
	local mrow=1
	forvalues i=5(1)15 {
	 quiet sum sim`x'90_`i'
	 matrix mcarloresults90[`mrow',`mcol']=r(mean) 
	 quiet sum sim`x'95_`i'
	 matrix mcarloresults95[`mrow',`mcol']=r(mean) 
	 quiet sum sim`x'99_`i'
	 matrix mcarloresults99[`mrow',`mcol']=r(mean) 
	 local mrow=`mrow'+1
	 }
	local mcol=`mcol'+1
 }
 }
sum simest03
sum simr203
gen avicc03= (simicc103+simicc203)/2
sum avicc03

sum simest04
sum simr204
gen avicc04= (simicc104+simicc204)/2
sum avicc04

sum simest05
sum simr205
gen avicc05= (simicc105+simicc205)/2
sum avicc05
 
matrix nobs=(50, 60, 70, 80, 90, 100, 110, 120, 130, 140, 150)'
matrix ncontests=(225, 270, 315, 360, 405, 450, 495, 540, 585, 630, 675)'
matrix line8=(.8,.8,.8,.8,.8,.8,.8,.8,.8,.8,.8)'

matrix A==J(11,3,1)
matrix power90=A -mcarloresults90
matrix power95=A -mcarloresults95
matrix power99=A -mcarloresults99
matrix highnoise=(nobs,ncontests, power90, power95 ,power99)
matrix colnames highnoise =  N Contests .3p90 .4p90 .5p90 .3p95 .4p95 .5p95 .3p99 .4p99 .5p99
matrix list highnoise

* power graph *
matrix x3l95=lownoise[1..11,6]
matrix x4l95=lownoise[1..11,7]
matrix x5l95=lownoise[1..11,8]

matrix x3h95=highnoise[1..11,6]
matrix x4h95=highnoise[1..11,7]
matrix x5h95=highnoise[1..11,8]

capture drop nobs1 
capture drop ncontests1 
capture drop x3l951 
capture drop x4l951 
capture drop x5l951 
capture drop x3h951 
capture drop x4h951 
capture drop x5h951
capture drop line8

local xlist "nobs ncontests x3l95 x4l95 x5l95 x3h95 x4h95 x5h95 line8"

foreach x of local xlist{
svmat double `x', name(`x')
}
}


cd "./results"

twoway  (line  x3l951 x4l951 x5l951 nobs1 , lpattern(longdash solid solid) lcolor(gs10 gs10 black) lwidth(thick thick thick) xaxis(1)) ///
 (line line8 ncontests1 , lpattern(blank)  xaxis(2)) ///
||,  xlabel(50(25)150, gmax angle(horizontal) axis(1)) ///
xlabel(225(75)675 , axis(2)) yscale(r(.5 1)) ylabel(.5(.1)1) ///
title("Low noise scenario") ytitle("Power") xtitle("Number of subjects", axis(1)) xtitle("Number of contests", axis(2)) legend(label(1 "-.3") label(2 "-.4") label(3 "-.5") order(3 2 1) pos(4) ring(0) region(color(none))) scheme(lean1) saving(powerl,replace) 

twoway  (line  x3h951 x4h951 x5h951 nobs1 , lpattern(longdash solid solid) lcolor(gs10 gs10 black) lwidth(thick thick thick) xaxis(1)) ///
 (line line8 ncontests1 , lpattern(blank)  xaxis(2)) ///
||,  xlabel(50(25)150, gmax angle(horizontal) axis(1)) ///
xlabel(225(75)675 , axis(2)) yscale(r(.5 1)) ylabel(.5(.1)1) ///
title("High noise scenario") ytitle("Power") xtitle("Number of subjects", axis(1)) xtitle("Number of contests", axis(2)) legend(label(1 "-.3") label(2 "-.4") label(3 "-.5") order(3 2 1) pos(4) ring(0) region(color(none))) scheme(lean1) saving(powerh,replace) 

gr combine powerl.gph powerh.gph, ///
col(2) iscale(0.7) xcommon scheme(lean1) 
graph export figs5.pdf, fontface(Arial) replace


**************************************************************************************
*Figure S6 Assumed data generating process for |Difference in strength|. 
**************************************************************************************
cd "..\"	
	 
use data_study4.dta, clear
set more off

cd "./results"
*abs.diff.in strength is folded normal
graph twoway (histogram absdstrength, width(0.1)  fcolor(gs10) ytitle("Density") xtitle("|Difference in strength|") xscale(titlegap(3)) yscale(titlegap(2)) scheme(lean1)) ///
(function y=(sqrt(2/_pi))*exp(-x^2), range(0 4) lpattern(solid) lcolor(black) lwidth(vthick)  ///
legend(label(1 "Actual data (Study 2)") label(2 "Simulation (folded normal distribution)")  pos(2) ring(0) region(color(none)) )), /// 
 saving(foldednormal, replace) 
graph export figs6.pdf, fontface(Arial) replace
 
**************************************************************************************
*Figure S7: Data generating process for contest duration with coefficient -.5
**************************************************************************************
cd "..\"	
*simulate 10,710 contests for scenario (low noise with coefficient -.5) 
quiet{

use startmontecarlo, clear
set more off
set matsize 800
set seed 1234567
*matrix montecarloDGP=J(765,10,0)
		*data generating process
		**mixture mean model
		gen a=2.5
		gen inva=1/a

forvalues k=1(1)14 {
		*generate absolute difference in strength according to a folded normal distribution (absolute value of difference of two standard normal random variables)
		by simid1, sort: gen u1=uniform()
		gen simstr1=invnormal(u1)
		by simid2, sort: gen u2=uniform()
		gen simstr2=invnormal(u2)
		capture drop simabsdstr
		gen simabsdstr=abs(simstr1-simstr2)
		
		*contestant random effects 
		quiet{
		capture drop a1 a2 
		gen a1=0
		gen a2=0


		forvalues i=1(1)170 {
			local aval=rnormal(0,1/2)
			replace a1=`aval' if simid1==`i' 
			replace a2=`aval' if simid2==`i' 
			}
		}
		
		
		** coefficient -.5 **
		gen xb05=3-.5*simabsdstr+a1+a2
		gen exb05=exp(xb05)
		*common error component for contest
		gen xg05=rgamma(inva,a)
		gen xbg05=exb05*xg05
		gen nby05=rpoisson(xbg05)
		
		*store in column vector
		mkmat nby05
		drop xb05 exb05 xg05 xbg05 nby05
		drop u1 u2 simstr1 simstr2

		if `k'==1 {
		matrix montecarloDGP=nby05
		}
		else {
		matrix montecarloDGP=(montecarloDGP,nby05)
		}
		
}
			svmat montecarloDGP
preserve
stack montecarloDGP*, into(durationsim) clear
gen participantid=0
save montecarloDGP	, replace		
restore

use data_study3,clear
set more off
merge m:m participantid using montecarloDGP.dta
gen cdurationsim=durationsim 
replace cdurationsim=225 if durationsim>=225
}

cd "./results"
histogram finalcosts, width(0.1)  fcolor(gs10) ytitle("Density") xtitle("Contest duration") xscale(titlegap(3)) yscale(titlegap(2)) scheme(lean1) title("Data from Study 2") saving(histact,replace)  
gen nsim=_n
histogram cdurationsim if nsim<10001, width(0.1)  fcolor(gs10) ytitle("Density") xtitle("Contest duration") xscale(titlegap(3)) yscale(titlegap(2)) scheme(lean1) title("Simulations (10,000 contests)") saving(histsim,replace) 

gr combine histact.gph histsim.gph, ///
col(2) iscale(0.7) xcommon scheme(lean1) 
graph export figs7.pdf, fontface(Arial) replace


*****************************************************************************************************************************
** RE-RUN Monte Carlo simulations **
** saves output to current directory (results)**
** WARNING: If the path "dir\matrices\" does not exist, matsave will result in the error **
*****************************************************************************************************************************

************************************************************************
***Monte Carlo preparation of panel structure ***	
************************************************************************
clear all
quiet{
* use seq function downloadable from dm44 from http://www.stata.com/stb/stb37
which seq
*matching matrix within session
set seed 1242342345
matrix a = (1,3,5,7,9,1,2,5,6,8,1,2,5,6,8,1,2,3,4,7,1,2,3,4,7,1,2,3,4,6,1,2,3,4,6,1,2,3,4,5,1,2,3,4,5\2,4,6,8,10,3,4,7,9,10,4,3,10,7,9,5,6,8,9,10,6,5,10,8,9,7,9,5,10,8,8,7,9,5,10,10,8,6,7,9,9,10,7,6,8)
matrix a=a'
matrix matching=a
set matsize 800
*stack matching matrix to duplicate for 16 sessions 
forvalues i=1(1)16{
matrix matching=matching\a
}

matrix drop a
set more off
svmat matching
matrix drop matching
*reset matsize to free up memory
set matsize 100

*"session number": generate sequence with blocks of size 45, starting at 1, ending at 22
seq simsess, f(1) t(17) b(45)
gen simid1=(simsess-1)*10+matching1
gen simid2=(simsess-1)*10+matching2
sort simsess 
drop matching*
gen simgameid=_n
}

save startmontecarlo, replace

************************************************************************
***Monte Carlo programs ***	
*modified version of Hilbe, 2011, Negative Binomial Regression, Cambridge, 2nd edition, Chapter 9.2, Table 9.6, p.229
*alpha=2.5, cons=4 
*run on beta coefficient sizes -.5, -.4, -.3
************************************************************************

**************************************************************
** Low noise variant **
**************************************************************

clear all
*Prepare a matrix montecarlo for results from Monte Carlo
matrix montecarlo03 = J(8,12,.)
matrix montecarlo04 = J(8,12,.)
matrix montecarlo05 = J(8,12,.)
*Columns: number of contests for number of sessions (with 10 participants each) 5 6 7 8 9 10 11 12 13 14 15 16
*matrix colnames montecarlo =  225 270 315 360 405 450 495 540 585 630 675 720
*matrix rownames montecarlo = estm pval reject90 reject95 reject99 r2 icc1 icc2
forvalues j=1(1)12 {
	forvalues i=1(1)8 {
	matrix montecarlo03[`i',`j']=0
	matrix montecarlo04[`i',`j']=0
	matrix montecarlo05[`i',`j']=0
	}
}


	capture program drop mymontecarlo
	program define mymontecarlo, rclass
		use startmontecarlo, clear
		*generate absolute difference in strength according to a folded normal distribution (absolute value of difference of two standard normal random variables)
		by simid1, sort: gen u1=uniform()
		gen simstr1=invnormal(u1)
		by simid2, sort: gen u2=uniform()
		gen simstr2=invnormal(u2)
		gen simabsdstr=abs(simstr1-simstr2)


		*contestant random effects 
		quiet{
		capture drop a1 a2 
		gen a1=0
		gen a2=0


		forvalues i=1(1)170 {
			local aval=rnormal(0,1/2)
			replace a1=`aval' if simid1==`i' 
			replace a2=`aval' if simid2==`i' 
			}
		}
		
		*data generating process
		**mixture mean model
		gen a=2.5
		gen inva=1/a
		** coefficient -.3 **
		gen xb03=3-.3*simabsdstr+a1+a2
		gen exb03=exp(xb03)
		*common error component for contest
		gen xg03=rgamma(inva,a)
		gen xbg03=exb03*xg03
		gen nby03=rpoisson(xbg03)
		
		** coefficient -.4 **
		gen xb04=3-.4*simabsdstr+a1+a2
		gen exb04=exp(xb04)
		*common error component for contest
		gen xg04=rgamma(inva,a)
		gen xbg04=exb04*xg04
		gen nby04=rpoisson(xbg04)
		
		** coefficient -.5 **
		gen xb05=3-.5*simabsdstr+a1+a2
		gen exb05=exp(xb05)
		*common error component for contest
		gen xg05=rgamma(inva,a)
		gen xbg05=exb05*xg05
		gen nby05=rpoisson(xbg05)
		
		

		preserve
			matsave montecarlo03, replace p("matrices") dropall
		restore
		preserve
			matsave montecarlo04, replace p("matrices") dropall
		restore
		preserve
			matsave montecarlo05, replace p("matrices") dropall
		restore
		


		** coefficient -.3 **
		local colno=1
		forvalues i=5(1)16 {
		matload montecarlo03,  p("matrices")  saving over
			preserve
				drop if simsess>`i'
				* OLS R2
				reg nby03 simabsdstr 
				matrix montecarlo03[6,`colno']=e(r2)
				vce2way nbreg nby03 simabsdstr , cluster(simid1 simid2)
				capture drop resid03 pred03
				predict double pred03, n
				gen resid03=nby03-pred03
				matrix temp=e(b)
				matrix montecarlo03[1,`colno']=temp[1,1]
				matrix drop temp
				test simabsdstr=0
				matrix temp=r(p)
				matrix montecarlo03[2,`colno']=temp[1,1]
				matrix drop temp
				* Rejection diff. strength coefficient==0
				scalar s90=0
				scalar s95=0
				scalar s99=0
				if r(p)>.1 {
					scalar s90=1
					scalar s95=1
					scalar s99=1
				}
				if r(p)>.05 {
					scalar s95=1
					scalar s99=1
				}
				if r(p)>.01 {
					scalar s99=1
				}
	
				matrix montecarlo03[3,`colno']=s90
				matrix montecarlo03[4,`colno']=s95
				matrix montecarlo03[5,`colno']=s99
				
				*Intracluster correlation coefficient
				loneway resid03 simid1
				matrix montecarlo03[7,`colno']=r(rho)
				loneway resid03 simid2
				matrix montecarlo03[8,`colno']=r(rho)
				
				
				matsave montecarlo03, replace p("matrices") dropall
			restore
			local colno=`colno'+1
			}

		
		** coefficient -.4 **
		local colno=1
		forvalues i=5(1)16 {
		matload montecarlo04,  p("matrices")  saving over
			preserve
				drop if simsess>`i'
				* OLS R2
				reg nby04 simabsdstr 
				matrix montecarlo04[6,`colno']=e(r2)
				vce2way nbreg nby04 simabsdstr , cluster(simid1 simid2)
				capture drop resid04 pred04
				predict double pred04, n
				gen resid04=nby04-pred04
				matrix temp=e(b)
				matrix montecarlo04[1,`colno']=temp[1,1]
				matrix drop temp
				test simabsdstr=0
				matrix temp=r(p)
				matrix montecarlo04[2,`colno']=temp[1,1]
				matrix drop temp
				* Rejection diff. strength coefficient==0
				scalar s90=0
				scalar s95=0
				scalar s99=0
				if r(p)>.1 {
					scalar s90=1
					scalar s95=1
					scalar s99=1
				}
				if r(p)>.05 {
					scalar s95=1
					scalar s99=1
				}
				if r(p)>.01 {
					scalar s99=1
				}
	
				matrix montecarlo04[3,`colno']=s90
				matrix montecarlo04[4,`colno']=s95
				matrix montecarlo04[5,`colno']=s99
				
				*Intracluster correlation coefficient
				loneway resid04 simid1
				matrix montecarlo04[7,`colno']=r(rho)
				loneway resid04 simid2
				matrix montecarlo04[8,`colno']=r(rho)
				
				
				matsave montecarlo04, replace p("matrices") dropall
			restore
			local colno=`colno'+1
			}

		** coefficient -.5 **
		local colno=1
		forvalues i=5(1)16 {
		matload montecarlo05,  p("matrices")  saving over
			preserve
				drop if simsess>`i'
				* OLS R2
				reg nby05 simabsdstr 
				matrix montecarlo05[6,`colno']=e(r2)
				vce2way nbreg nby05 simabsdstr , cluster(simid1 simid2)
				capture drop resid05 pred05
				predict double pred05, n
				gen resid05=nby05-pred05
				matrix temp=e(b)
				matrix montecarlo05[1,`colno']=temp[1,1]
				matrix drop temp
				test simabsdstr=0
				matrix temp=r(p)
				matrix montecarlo05[2,`colno']=temp[1,1]
				matrix drop temp
				* Rejection diff. strength coefficient==0
				scalar s90=0
				scalar s95=0
				scalar s99=0
				if r(p)>.1 {
					scalar s90=1
					scalar s95=1
					scalar s99=1
				}
				if r(p)>.05 {
					scalar s95=1
					scalar s99=1
				}
				if r(p)>.01 {
					scalar s99=1
				}
	
				matrix montecarlo05[3,`colno']=s90
				matrix montecarlo05[4,`colno']=s95
				matrix montecarlo05[5,`colno']=s99
				
				*Intracluster correlation coefficient
				loneway resid05 simid1
				matrix montecarlo05[7,`colno']=r(rho)
				loneway resid05 simid2
				matrix montecarlo05[8,`colno']=r(rho)
				
				
				matsave montecarlo05, replace p("matrices") dropall
			restore
			local colno=`colno'+1
			}
			
			
			
			
			
			
			
			
			
			*return values	
		local colno=1
		forvalues i=5(1)15 {
			return scalar sim0390_`i'=montecarlo03[3,`colno']
			return scalar sim0395_`i'=montecarlo03[4,`colno']
			return scalar sim0399_`i'=montecarlo03[5,`colno']
			return scalar sim0490_`i'=montecarlo04[3,`colno']
			return scalar sim0495_`i'=montecarlo04[4,`colno']
			return scalar sim0499_`i'=montecarlo04[5,`colno']
			return scalar sim0590_`i'=montecarlo05[3,`colno']
			return scalar sim0595_`i'=montecarlo05[4,`colno']
			return scalar sim0599_`i'=montecarlo05[5,`colno']
			local colno=`colno'+1
			}
			return scalar simest03=montecarlo03[1,8]
			return scalar simr203=montecarlo03[6,8]
			return scalar simicc103=montecarlo03[7,8]
			return scalar simicc203=montecarlo03[8,8]
			
			return scalar simest04=montecarlo04[1,8]
			return scalar simr204=montecarlo04[6,8]
			return scalar simicc104=montecarlo04[7,8]
			return scalar simicc204=montecarlo04[8,8]
			
			return scalar simest05=montecarlo05[1,8]
			return scalar simr205=montecarlo05[6,8]
			return scalar simicc105=montecarlo05[7,8]
			return scalar simicc205=montecarlo05[8,8]
			

	end		


	*Run Monte Carlo
	set more off
	preserve
		simulate simest03=r(simest03) simr203=r(simr203) simicc103=r(simicc103) simicc203=r(simicc203) ///
		sim0390_5=r(sim0390_5) sim0395_5=r(sim0395_5) sim0399_5=r(sim0399_5) sim0390_6=r(sim0390_6) sim0395_6=r(sim0395_6) sim0399_6=r(sim0399_6) ///
		sim0390_7=r(sim0390_7) sim0395_7=r(sim0395_7) sim0399_7=r(sim0399_7) sim0390_8=r(sim0390_8) sim0395_8=r(sim0395_8) sim0399_8=r(sim0399_8) ///
		sim0390_9=r(sim0390_9) sim0395_9=r(sim0395_9) sim0399_9=r(sim0399_9) sim0390_10=r(sim0390_10) sim0395_10=r(sim0395_10) sim0399_10=r(sim0399_10) ///
		sim0390_11=r(sim0390_11) sim0395_11=r(sim0395_11) sim0399_11=r(sim0399_11) sim0390_12=r(sim0390_12) sim0395_12=r(sim0395_12) sim0399_12=r(sim0399_12) ///
		sim0390_13=r(sim0390_13) sim0395_13=r(sim0395_13) sim0399_13=r(sim0399_13) sim0390_14=r(sim0390_14) sim0395_14=r(sim0395_14) sim0399_14=r(sim0399_14) ///
		sim0390_15=r(sim0390_15) sim0395_15=r(sim0395_15) sim0399_15=r(sim0399_15)  ///
		simest04=r(simest04) simr204=r(simr204) simicc104=r(simicc104) simicc204=r(simicc204) ///
		sim0490_5=r(sim0490_5) sim0495_5=r(sim0495_5) sim0499_5=r(sim0499_5) sim0490_6=r(sim0490_6) sim0495_6=r(sim0495_6) sim0499_6=r(sim0499_6) ///
		sim0490_7=r(sim0490_7) sim0495_7=r(sim0495_7) sim0499_7=r(sim0499_7) sim0490_8=r(sim0490_8) sim0495_8=r(sim0495_8) sim0499_8=r(sim0499_8) ///
		sim0490_9=r(sim0490_9) sim0495_9=r(sim0495_9) sim0499_9=r(sim0499_9) sim0490_10=r(sim0490_10) sim0495_10=r(sim0495_10) sim0499_10=r(sim0499_10) ///
		sim0490_11=r(sim0490_11) sim0495_11=r(sim0495_11) sim0499_11=r(sim0499_11) sim0490_12=r(sim0490_12) sim0495_12=r(sim0495_12) sim0499_12=r(sim0499_12) ///
		sim0490_13=r(sim0490_13) sim0495_13=r(sim0495_13) sim0499_13=r(sim0499_13) sim0490_14=r(sim0490_14) sim0495_14=r(sim0495_14) sim0499_14=r(sim0499_14) ///
		sim0490_15=r(sim0490_15) sim0495_15=r(sim0495_15) sim0499_15=r(sim0499_15)  ///
		simest05=r(simest05) simr205=r(simr205) simicc105=r(simicc105) simicc205=r(simicc205) ///
		sim0590_5=r(sim0590_5) sim0595_5=r(sim0595_5) sim0599_5=r(sim0599_5) sim0590_6=r(sim0590_6) sim0595_6=r(sim0595_6) sim0599_6=r(sim0599_6) ///
		sim0590_7=r(sim0590_7) sim0595_7=r(sim0595_7) sim0599_7=r(sim0599_7) sim0590_8=r(sim0590_8) sim0595_8=r(sim0595_8) sim0599_8=r(sim0599_8) ///
		sim0590_9=r(sim0590_9) sim0595_9=r(sim0595_9) sim0599_9=r(sim0599_9) sim0590_10=r(sim0590_10) sim0595_10=r(sim0595_10) sim0599_10=r(sim0599_10) ///
		sim0590_11=r(sim0590_11) sim0595_11=r(sim0595_11) sim0599_11=r(sim0599_11) sim0590_12=r(sim0590_12) sim0595_12=r(sim0595_12) sim0599_12=r(sim0599_12) ///
		sim0590_13=r(sim0590_13) sim0595_13=r(sim0595_13) sim0599_13=r(sim0599_13) sim0590_14=r(sim0590_14) sim0595_14=r(sim0595_14) sim0599_14=r(sim0599_14) ///
		sim0590_15=r(sim0590_15) sim0595_15=r(sim0595_15) sim0599_15=r(sim0599_15) , ///
		reps(1000) seed(12345) saving(mymontecarlolownoise, replace): mymontecarlo
	 restore

	 
**************************************************************
** high noise variant **
**************************************************************

clear all
*Prepare a matrix montecarlo for results from Monte Carlo
matrix montecarlo03 = J(8,12,.)
matrix montecarlo04 = J(8,12,.)
matrix montecarlo05 = J(8,12,.)
*Columns: number of contests for number of sessions (with 10 participants each) 5 6 7 8 9 10 11 12 13 14 15 16
*matrix colnames montecarlo =  225 270 315 360 405 450 495 540 585 630 675 720
*matrix rownames montecarlo = estm pval reject90 reject95 reject99 r2 icc1 icc2
forvalues j=1(1)12 {
	forvalues i=1(1)8 {
	matrix montecarlo03[`i',`j']=0
	matrix montecarlo04[`i',`j']=0
	matrix montecarlo05[`i',`j']=0
	}
}


	capture program drop mymontecarlo
	program define mymontecarlo, rclass
		use startmontecarlo, clear
		*generate absolute difference in strength according to a folded normal distribution (absolute value of difference of two standard normal random variables)
		by simid1, sort: gen u1=uniform()
		gen simstr1=invnormal(u1)
		by simid2, sort: gen u2=uniform()
		gen simstr2=invnormal(u2)
		gen simabsdstr=abs(simstr1-simstr2)


		*contestant random effects 
		quiet{
		capture drop a1 a2 
		gen a1=0
		gen a2=0


		forvalues i=1(1)170 {
			local aval=rnormal(0,3/4)
			replace a1=`aval' if simid1==`i' 
			replace a2=`aval' if simid2==`i' 
			}
		}
		
		*data generating process
		**mixture mean model
		gen a=2.5
		gen inva=1/a
		** coefficient -.3 **
		gen xb03=3-.3*simabsdstr+a1+a2
		gen exb03=exp(xb03)
		*common error component for contest
		gen xg03=rgamma(inva,a)
		gen xbg03=exb03*xg03
		gen nby03=rpoisson(xbg03)
		
		** coefficient -.4 **
		gen xb04=3-.4*simabsdstr+a1+a2
		gen exb04=exp(xb04)
		*common error component for contest
		gen xg04=rgamma(inva,a)
		gen xbg04=exb04*xg04
		gen nby04=rpoisson(xbg04)
		
		** coefficient -.5 **
		gen xb05=3-.5*simabsdstr+a1+a2
		gen exb05=exp(xb05)
		*common error component for contest
		gen xg05=rgamma(inva,a)
		gen xbg05=exb05*xg05
		gen nby05=rpoisson(xbg05)
		
		

		preserve
			matsave montecarlo03, replace p("matrices") dropall
		restore
		preserve
			matsave montecarlo04, replace p("matrices") dropall
		restore
		preserve
			matsave montecarlo05, replace p("matrices") dropall
		restore
		


		** coefficient -.3 **
		local colno=1
		forvalues i=5(1)16 {
		matload montecarlo03,  p("matrices")  saving over
			preserve
				drop if simsess>`i'
				* OLS R2
				reg nby03 simabsdstr 
				matrix montecarlo03[6,`colno']=e(r2)
				vce2way nbreg nby03 simabsdstr , cluster(simid1 simid2)
				capture drop resid03 pred03
				predict double pred03, n
				gen resid03=nby03-pred03
				matrix temp=e(b)
				matrix montecarlo03[1,`colno']=temp[1,1]
				matrix drop temp
				test simabsdstr=0
				matrix temp=r(p)
				matrix montecarlo03[2,`colno']=temp[1,1]
				matrix drop temp
				* Rejection diff. strength coefficient==0
				scalar s90=0
				scalar s95=0
				scalar s99=0
				if r(p)>.1 {
					scalar s90=1
					scalar s95=1
					scalar s99=1
				}
				if r(p)>.05 {
					scalar s95=1
					scalar s99=1
				}
				if r(p)>.01 {
					scalar s99=1
				}
	
				matrix montecarlo03[3,`colno']=s90
				matrix montecarlo03[4,`colno']=s95
				matrix montecarlo03[5,`colno']=s99
				
				*Intracluster correlation coefficient
				loneway resid03 simid1
				matrix montecarlo03[7,`colno']=r(rho)
				loneway resid03 simid2
				matrix montecarlo03[8,`colno']=r(rho)
				
				
				matsave montecarlo03, replace p("matrices") dropall
			restore
			local colno=`colno'+1
			}

		
		** coefficient -.4 **
		local colno=1
		forvalues i=5(1)16 {
		matload montecarlo04,  p("matrices")  saving over
			preserve
				drop if simsess>`i'
				* OLS R2
				reg nby04 simabsdstr 
				matrix montecarlo04[6,`colno']=e(r2)
				vce2way nbreg nby04 simabsdstr , cluster(simid1 simid2)
				capture drop resid04 pred04
				predict double pred04, n
				gen resid04=nby04-pred04
				matrix temp=e(b)
				matrix montecarlo04[1,`colno']=temp[1,1]
				matrix drop temp
				test simabsdstr=0
				matrix temp=r(p)
				matrix montecarlo04[2,`colno']=temp[1,1]
				matrix drop temp
				* Rejection diff. strength coefficient==0
				scalar s90=0
				scalar s95=0
				scalar s99=0
				if r(p)>.1 {
					scalar s90=1
					scalar s95=1
					scalar s99=1
				}
				if r(p)>.05 {
					scalar s95=1
					scalar s99=1
				}
				if r(p)>.01 {
					scalar s99=1
				}
	
				matrix montecarlo04[3,`colno']=s90
				matrix montecarlo04[4,`colno']=s95
				matrix montecarlo04[5,`colno']=s99
				
				*Intracluster correlation coefficient
				loneway resid04 simid1
				matrix montecarlo04[7,`colno']=r(rho)
				loneway resid04 simid2
				matrix montecarlo04[8,`colno']=r(rho)
				
				
				matsave montecarlo04, replace p("matrices") dropall
			restore
			local colno=`colno'+1
			}

		** coefficient -.5 **
		local colno=1
		forvalues i=5(1)16 {
		matload montecarlo05,  p("matrices")  saving over
			preserve
				drop if simsess>`i'
				* OLS R2
				reg nby05 simabsdstr 
				matrix montecarlo05[6,`colno']=e(r2)
				vce2way nbreg nby05 simabsdstr , cluster(simid1 simid2)
				capture drop resid05 pred05
				predict double pred05, n
				gen resid05=nby05-pred05
				matrix temp=e(b)
				matrix montecarlo05[1,`colno']=temp[1,1]
				matrix drop temp
				test simabsdstr=0
				matrix temp=r(p)
				matrix montecarlo05[2,`colno']=temp[1,1]
				matrix drop temp
				* Rejection diff. strength coefficient==0
				scalar s90=0
				scalar s95=0
				scalar s99=0
				if r(p)>.1 {
					scalar s90=1
					scalar s95=1
					scalar s99=1
				}
				if r(p)>.05 {
					scalar s95=1
					scalar s99=1
				}
				if r(p)>.01 {
					scalar s99=1
				}
	
				matrix montecarlo05[3,`colno']=s90
				matrix montecarlo05[4,`colno']=s95
				matrix montecarlo05[5,`colno']=s99
				
				*Intracluster correlation coefficient
				loneway resid05 simid1
				matrix montecarlo05[7,`colno']=r(rho)
				loneway resid05 simid2
				matrix montecarlo05[8,`colno']=r(rho)
				
				
				matsave montecarlo05, replace p("matrices") dropall
			restore
			local colno=`colno'+1
			}
			
			
			
			
			
			
			
			
			
			*return values	
		local colno=1
		forvalues i=5(1)15 {
			return scalar sim0390_`i'=montecarlo03[3,`colno']
			return scalar sim0395_`i'=montecarlo03[4,`colno']
			return scalar sim0399_`i'=montecarlo03[5,`colno']
			return scalar sim0490_`i'=montecarlo04[3,`colno']
			return scalar sim0495_`i'=montecarlo04[4,`colno']
			return scalar sim0499_`i'=montecarlo04[5,`colno']
			return scalar sim0590_`i'=montecarlo05[3,`colno']
			return scalar sim0595_`i'=montecarlo05[4,`colno']
			return scalar sim0599_`i'=montecarlo05[5,`colno']
			local colno=`colno'+1
			}
			return scalar simest03=montecarlo03[1,8]
			return scalar simr203=montecarlo03[6,8]
			return scalar simicc103=montecarlo03[7,8]
			return scalar simicc203=montecarlo03[8,8]
			
			return scalar simest04=montecarlo04[1,8]
			return scalar simr204=montecarlo04[6,8]
			return scalar simicc104=montecarlo04[7,8]
			return scalar simicc204=montecarlo04[8,8]
			
			return scalar simest05=montecarlo05[1,8]
			return scalar simr205=montecarlo05[6,8]
			return scalar simicc105=montecarlo05[7,8]
			return scalar simicc205=montecarlo05[8,8]
			

	end		


	*Run Monte Carlo
	set more off
	preserve
		simulate simest03=r(simest03) simr203=r(simr203) simicc103=r(simicc103) simicc203=r(simicc203) ///
		sim0390_5=r(sim0390_5) sim0395_5=r(sim0395_5) sim0399_5=r(sim0399_5) sim0390_6=r(sim0390_6) sim0395_6=r(sim0395_6) sim0399_6=r(sim0399_6) ///
		sim0390_7=r(sim0390_7) sim0395_7=r(sim0395_7) sim0399_7=r(sim0399_7) sim0390_8=r(sim0390_8) sim0395_8=r(sim0395_8) sim0399_8=r(sim0399_8) ///
		sim0390_9=r(sim0390_9) sim0395_9=r(sim0395_9) sim0399_9=r(sim0399_9) sim0390_10=r(sim0390_10) sim0395_10=r(sim0395_10) sim0399_10=r(sim0399_10) ///
		sim0390_11=r(sim0390_11) sim0395_11=r(sim0395_11) sim0399_11=r(sim0399_11) sim0390_12=r(sim0390_12) sim0395_12=r(sim0395_12) sim0399_12=r(sim0399_12) ///
		sim0390_13=r(sim0390_13) sim0395_13=r(sim0395_13) sim0399_13=r(sim0399_13) sim0390_14=r(sim0390_14) sim0395_14=r(sim0395_14) sim0399_14=r(sim0399_14) ///
		sim0390_15=r(sim0390_15) sim0395_15=r(sim0395_15) sim0399_15=r(sim0399_15)  ///
		simest04=r(simest04) simr204=r(simr204) simicc104=r(simicc104) simicc204=r(simicc204) ///
		sim0490_5=r(sim0490_5) sim0495_5=r(sim0495_5) sim0499_5=r(sim0499_5) sim0490_6=r(sim0490_6) sim0495_6=r(sim0495_6) sim0499_6=r(sim0499_6) ///
		sim0490_7=r(sim0490_7) sim0495_7=r(sim0495_7) sim0499_7=r(sim0499_7) sim0490_8=r(sim0490_8) sim0495_8=r(sim0495_8) sim0499_8=r(sim0499_8) ///
		sim0490_9=r(sim0490_9) sim0495_9=r(sim0495_9) sim0499_9=r(sim0499_9) sim0490_10=r(sim0490_10) sim0495_10=r(sim0495_10) sim0499_10=r(sim0499_10) ///
		sim0490_11=r(sim0490_11) sim0495_11=r(sim0495_11) sim0499_11=r(sim0499_11) sim0490_12=r(sim0490_12) sim0495_12=r(sim0495_12) sim0499_12=r(sim0499_12) ///
		sim0490_13=r(sim0490_13) sim0495_13=r(sim0495_13) sim0499_13=r(sim0499_13) sim0490_14=r(sim0490_14) sim0495_14=r(sim0495_14) sim0499_14=r(sim0499_14) ///
		sim0490_15=r(sim0490_15) sim0495_15=r(sim0495_15) sim0499_15=r(sim0499_15)  ///
		simest05=r(simest05) simr205=r(simr205) simicc105=r(simicc105) simicc205=r(simicc205) ///
		sim0590_5=r(sim0590_5) sim0595_5=r(sim0595_5) sim0599_5=r(sim0599_5) sim0590_6=r(sim0590_6) sim0595_6=r(sim0595_6) sim0599_6=r(sim0599_6) ///
		sim0590_7=r(sim0590_7) sim0595_7=r(sim0595_7) sim0599_7=r(sim0599_7) sim0590_8=r(sim0590_8) sim0595_8=r(sim0595_8) sim0599_8=r(sim0599_8) ///
		sim0590_9=r(sim0590_9) sim0595_9=r(sim0595_9) sim0599_9=r(sim0599_9) sim0590_10=r(sim0590_10) sim0595_10=r(sim0595_10) sim0599_10=r(sim0599_10) ///
		sim0590_11=r(sim0590_11) sim0595_11=r(sim0595_11) sim0599_11=r(sim0599_11) sim0590_12=r(sim0590_12) sim0595_12=r(sim0595_12) sim0599_12=r(sim0599_12) ///
		sim0590_13=r(sim0590_13) sim0595_13=r(sim0595_13) sim0599_13=r(sim0599_13) sim0590_14=r(sim0590_14) sim0595_14=r(sim0595_14) sim0599_14=r(sim0599_14) ///
		sim0590_15=r(sim0590_15) sim0595_15=r(sim0595_15) sim0599_15=r(sim0599_15) , ///
		reps(1000) seed(12345) saving(mymontecarlohighnoise, replace): mymontecarlo
	 restore

******************************************************************************************************************************************************************************
******************************************************************************************************************************************************************************


	
**********************************************************
***********Study 3 (uses saved results from Study 2)***********
**********************************************************
cd "..\"	
	 
use data_study3.dta, clear
set more off

cd "./results"

**********************************************************
***********Prepare estimation ***********
**********************************************************
	

label var absdstrength "|Difference in strength|"
label var absdheight "|Difference in height|"
label var absdstrength2 "|Difference in strength|^2"


*figure scheme
set scheme lean1

*Note: Because one subject was in game with excluded subject (shoulder injury) in every period, can't use fixed period to get stats for all
by participantid, sort: egen uniqueperiod=min(period)

********************************************************************************
*********** Stats/results in manuscript text ***********
********************************************************************************


sum age if period==uniqueperiod

** Cronbach Alpha for composite strength **
alpha z_flexedbicep z_maxgrip z_maxchest z_perceivedstrength if period==uniqueperiod, item

** Conflict duration **
sum finalcosts if treatment==1 & contestantincluded==1 //MAT


** Regression coefficients reported / marginal effects **
*DURATION
quiet: sum finalcosts if treatment==1 & contestantincluded==1
local sdhat1=r(sd)
quiet: sum absdstrength if contestantincluded==1
local sdhat2=r(sd)
preserve 
keep if treatment==1 & contestantincluded==1 
vce2way nbreg finalcosts absdstrength  per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 per11 per12 per13, cluster(id1 id2)
** standardized average marginal effect **
margins, dydx(*) post
matrix temp=e(b)
matrix list temp
local bhat=temp[1,1]
local effectsize=(`bhat'*`sdhat2')/`sdhat1'
display `effectsize' "(standardized average marginal effect)"
restore

*WINNER
*Test that absdstrength and absdstrength2 are jointly 0 in regression for winner*
preserve
keep if treatment==1 & contestantincluded==1 & draw!=1
quietly: vce2way logit strongwins absdstrength c.absdstrength#c.absdstrength per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 per11 per12 per13 , cluster(id1 id2)
test absdstrength c.absdstrength#c.absdstrength
restore
*Discrete marginal effect on prob(strongest wins)
*percentage point increase in the probability of the stronger contestant winning when the difference in strength increases by one standard deviation from its mean
preserve
quiet: sum absdstrength if contestantincluded==1
local sdhat2=r(sd)
keep if treatment==1 & contestantincluded==1 
vce2way logit strongwins absdstrength c.absdstrength#c.absdstrength  per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 per11 per12 per13, cluster(id1 id2)
mtable, atmeans at(absdstrength = gen(absdstrength))  ///
		at(absdstrength = gen(absdstrength+`sdhat2'))  post
mlincom 2 - 1, rowname(1std.dev.increase |diff.strength|)
matrix temp=_mlincom
matrix list temp
local effectsize=temp[1,1]
display `effectsize' "(average change in probability that stronger contestant wins for one std.dev. change in |difference strength|)"
restore
*Discrete marginal effect on prob(strongest wins), zero duration contests
preserve
quiet: sum absdstrength if contestantincluded==1
local mean2=r(mean)
local sdhat2=r(sd)
local mean2plus= `mean2'+`sdhat2'
keep if treatment==1 & contestantincluded==1 & finalcosts==0
vce2way logit strongwins absdstrength per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 per11 per12 per13, cluster(id1 id2)
mtable, atmeans at(absdstrength = gen(absdstrength))  ///
		at(absdstrength = gen(absdstrength+`sdhat2'))  post
mlincom 2 - 1, rowname(1std.dev.increase |diff.strength|)
matrix temp=_mlincom
matrix list temp
local effectsize=temp[1,1]
display `effectsize' "(average change in probability that stronger contestant wins for one std.dev. change in |difference strength|; zero duration contests)"
restore





********************************************************************************
***********Fig. 2 (complete): The effect of differences in strength on conflict duration in the Mutual Assessment Treatment***********
********************************************************************************
*lines in plot area for 5th and 95th percentile of absdstrength
sum absdstrength
centile absdstrength, centile(5 95)

preserve 
keep if treatment==1 & contestantincluded==1 
quietly: vce2way nbreg finalcosts absdstrength  per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 per11 per12 per13, cluster(id1 id2)
quietly: margins, at(absdstrength=(0(0.1)3.9)) atmeans
marginsplot,xline(.0639445) xline(1.889784) recast(line) recastci(rline) plotopt(lw(medthick) lcolor(dknavy)) ciopt(lpat(dash) lw(vthin)) ytitle("Predicted duration in seconds", size(small)) xtitle("Absolute difference in strength", size(small))  title("", size(small))  xscale(titlegap(3))  ylabel(0(20)80) yscale(r(-5,90) titlegap(2)) scheme(lean1) saving(durationSTRonly4, replace) 
restore

*Combine figures from studies 2 & 3
gr combine durationSTRonly.gph durationSTRonly4.gph, ///
	col(2) iscale(.7) xcommon scheme(lean1) altshrink  title("Study 2                                                                                   Study 3", size(small)) 
graph export figure2.pdf, fontface(Arial) replace



********************************************************************************
***********Fig. 3 (complete): Who wins the contests in the Mutual Assessment Treatment***********
********************************************************************************
*All contests
preserve
keep if treatment==1 & contestantincluded==1 & draw!=1
quietly: vce2way logit strongwins absdstrength c.absdstrength#c.absdstrength per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 per11 per12 per13 , cluster(id1 id2)
quietly: margins, at(absdstrength=(0(0.1)3.9)) atmeans
marginsplot,xline(.0639445) xline(1.889784) recast(line) recastci(rline) plotopt(lw(medthick) lcolor(dknavy)) ciopt(lpat(dash) lw(vthin)) ytitle("Probability of stronger winning", size(small)) xtitle("Absolute difference in strength", size(small)) title("b) Strongest wins", size(small))  xscale(titlegap(3)) ylabel(0(.2)1) yscale(r(-.3,1) titlegap(2))  scheme(lean1) saving(winner4,replace)  
restore
*zero duration contests
preserve
keep if treatment==1 & contestantincluded==1 & finalcosts==0
quietly: vce2way logit strongwins absdstrength per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 per11 per12 per13, cluster(id1 id2)
quietly: margins, at(absdstrength=(0(0.1)3.9)) atmeans
marginsplot,xline(.0639445) xline(1.889784) recast(line) recastci(rline) plotopt(lw(medthick) lcolor(dknavy)) ciopt(lpat(dash) lw(vthin)) ytitle("Probability of stronger winning", size(small)) xtitle("Absolute difference in strength", size(small)) title("d) Strongest wins (zero duration contests)", size(small))  xscale(titlegap(3)) ylabel(0(.2)1) yscale(r(-.3,1.2) titlegap(2)) scheme(lean1) saving(winnerzero4,replace)  
restore

*Fig 3
*Combine figures from studies 2 & 4
gr combine  winner.gph winner4.gph winnerzero.gph winnerzero4.gph, ///
	col(2) iscale(.7) xcommon scheme(lean1) altshrink  title("Study 2                                                                                   Study 3", size(small)) 
graph export figure3.pdf, fontface(Arial) replace







***********************************************
***********Fig. 4 (complete): Evidence for mutual assessment***********
***********************************************
gen dummyregressor1=finalcosts
gen dummyregressor2=-finalcosts


*Joint test of both coefficients being zero
preserve 
keep if treatment==1 & contestantincluded==1 
vce2way nbreg finalcosts strongformid weakformid per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 , cluster(id1 id2)
test strongformid_z  weakformid_z
restore

preserve 
keep if treatment==1 & contestantincluded==1 
quiet: vce2way nbreg finalcosts strongformid weakformid per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 per11 per12 per13 , cluster(id1 id2)
quiet: margins, at(weakformid=(-1.61(0.1)0.831)) atmeans
quiet: marginsplot, recast(line) recastci(rline) ciopts(lpattern(dash)) plotopt(lcolor(dknavy) lw(thick)) xlabel(#5) legend(off) scheme(lean1) ytitle("Predicted duration in seconds") xtitle("Strength of weaker contestant" "(standardized)") title("c) Observed results in Study 3" " ") xscale(titlegap(3)) yscale(range(0 150) titlegap(2)) ylabel(0(50)150) saving(weakformid4,replace)
quiet: marginsplot, recast(line) recastci(rline) ciopts(lpattern(dash)) plotopt(lcolor(dknavy) lw(thick)) xlabel(#5) legend(off) scheme(lean1) ytitle("Predicted duration in seconds") xtitle("Strength of weaker contestant" "(standardized)") title("d) Observed results in Study 3" " ") xscale(titlegap(3)) yscale(range(0 150) titlegap(2)) ylabel(0(50)150) saving(weakformid4d,replace)
quiet: vce2way nbreg finalcosts strongformid weakformid per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 per11 per12 per13 , cluster(id1 id2)
quiet: margins, at(strongformid=(-1.26(0.1)2.30)) atmeans
quietly: marginsplot, recast(line) recastci(rline) ciopts(lpattern(dash)) plotopt(lcolor(dknavy) lw(thick)) xlabel(#5) legend(off) scheme(lean1) ytitle("Predicted duration in seconds") xtitle("Strength of stronger contestant" "(standardized)") title("") yscale(range(0 150) titlegap(2)) ylabel(0(50)150) saving(strongformid4,replace)
graph twoway (function x/2, range(strongform) n(2) lpattern(dash)  lpattern(dash) lcolor(dknavy) lw(vthick) xlabel(#5) legend(off) scheme(lean1) ytitle("Duration") yscale(range(-3 3))  ylabel(none ) xlabel(none -1.2"Low" 2.2"High") title("a) Theoretical predictions:" "Mutual Assessment Model") xtitle("Strength of weaker contestant" "") saving(weak_dummy1,replace))
graph twoway (function -x/2, range(strongform) n(2) lpattern(dash)  lpattern(dash) lcolor(dknavy) lw(vthick) xlabel(#5) legend(off) scheme(lean1) ytitle("Duration") yscale(range(-3 3))  ylabel(none ) xlabel(none -1.2"Low" 2.2"High") title("") xtitle("Strength of stronger contestant" "") saving(strong_dummy1,replace))
graph twoway (function x/2, range(strongform) n(2) lpattern(dash)  lpattern(dash) lcolor(dknavy) lw(vthick) xlabel(#5) legend(off) scheme(lean1) ytitle("Duration") yscale(range(-3 3))  ylabel(none ) xlabel(none -1.2"Low" 2.2"High") title("b) Theoretical predictions:" "Pure Self-Assessment Model") xtitle("Strength of weaker contestant" "") saving(weak_dummy2,replace))
graph twoway (function x/4, range(strongform) n(2) lpattern(dash)  lpattern(dash) lcolor(dknavy) lw(vthick) xlabel(#5) legend(off) scheme(lean1) ytitle("Duration") yscale(range(-3 3))  ylabel(none ) xlabel(none -1.2"Low" 2.2"High") title("") xtitle("Strength of stronger contestant" "") saving(strong_dummy2,replace))
restore

*Combine figures from studies 2 & 3
gr combine weak_dummy1.gph weak_dummy2.gph weakformid.gph weakformid4d.gph strong_dummy1.gph strong_dummy2.gph strongformid.gph strongformid4.gph, col(4) iscale(0.5) scheme(lean1)  altshrink 
graph export figure4.pdf, fontface(Arial) replace

***************************************************************
***********Figure S4 Study 3: Distribution of the marginal effects ***********
***************************************************************
*distribution of marginal effects
sum finalcost if treatment==1 & contestantincluded==1
local sdhat1=r(sd)
sum absdstrength if contestantincluded==1
local sdhat2=r(sd)
preserve 
keep if treatment==1 & contestantincluded==1 
vce2way nbreg finalcosts absdstrength  per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 per11 per12 per13, cluster(id1 id2)
*standardized marginal effect 
margins, generate(MEduration) ///
at(absdstrength = gen(absdstrength)) at(absdstrength = gen(absdstrength+`sdhat2'))
gen DCduration = (MEduration2 - MEduration1)/`sdhat1'
hist DCduration, ytitle("Density") xtitle("Standardized marginal effect on contest duration") title("a) Contest duration") xscale(titlegap(3)) yscale(titlegap(2)) scheme(lean1) saving(MEduration, replace) 
sum DCduration
restore

*distribution of increase in prob(stronger wins) for one std.dev. increase in |diff strength|
preserve
quiet: sum absdstrength if contestantincluded==1
local sdhat2=r(sd)
keep if treatment==1 & contestantincluded==1 & draw!=1
vce2way logit strongwins absdstrength  c.absdstrength#c.absdstrength per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 per11 per12 per13, cluster(id1 id2)
*standardized marginal effect at means of covariates
margins, generate(MEwinner) ///
at(absdstrength = gen(absdstrength)) at(absdstrength = gen(absdstrength+`sdhat2'))
gen DCwinner = MEwinner2 - MEwinner1
hist DCwinner, ytitle("Density") xtitle("Change in prob(strongest wins)") title("b) Strongest wins") xscale(titlegap(3)) yscale(titlegap(2)) scheme(lean1) saving(MEwinner, replace) 
restore

*distribution of increase in prob(stronger wins) for one std.dev. increase in |diff strength|; for zero duration contests
preserve
quiet: sum absdstrength if contestantincluded==1
local sdhat2=r(sd)
keep if treatment==1 & contestantincluded==1 & finalcost==0
vce2way logit strongwins absdstrength per1 per2 per3 per4 per5 per6 per7 per8 per9 per10, cluster(id1 id2)
*standardized marginal effect at means of covariates
margins, generate(MEwinnerzero) ///
at(absdstrength = gen(absdstrength)) at(absdstrength = gen(absdstrength+`sdhat2'))
gen DCwinnerzero = MEwinnerzero2 - MEwinnerzero1
hist DCwinnerzero, ytitle("Density") xtitle("Change in prob(strongest wins)") title("c) Strongest wins (zero duration contests)") xscale(titlegap(3)) yscale(titlegap(2)) scheme(lean1) saving(MEwinnerzero, replace) 
restore

*Distribution of the effects of a one standard deviation" "increase in the absolute difference in strength
gr combine MEduration.gph MEwinner.gph MEwinnerzero.gph, ///
	col(2) hole(2) iscale(0.7) xcommon title("") scheme(lean1) 
graph export figureS4.pdf, fontface(Arial) replace







**********************************************************
***********Results (ESM) ***********
**********************************************************

*********** Number of contests (ESM) ***********
* tabulate the maximum number of periods against sessions: number of participants=number of periods+1
by sessionid period, sort: gen sessnum=_N
capture, drop sessnum
tab sessnum if period==uniqueperiod & treatment==1  //MAT

*********** Coefficients on the stronger contestant and weaker contestants’ strength (ESM) ***********
*the coefficients on the stronger contestant and weaker contestants’ strength have the same absolute value*

preserve 
*model 7
keep if treatment==1 & contestantincluded==1
vce2way nbreg finalcosts strongformid weakformid per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 per11 per12 per13  , cluster(id1 id2)
test strongformid+weakformid=0
restore


***********Correlation of Formidability and risk tolerance (ESM) ***********

*"we find no correlation between strength and risk tolerance"
preserve
collapse Formidability_z maxgrip maxchest flexedbicep perceivedstrength age height weight bmi  risk treatment, by(participantid)
pwcorr  risk Formidability_z, sig obs
restore 

**********************************************************
***********Table S1: Descriptive statistics for Study 3 (non-standardized) ***********
**********************************************************

*Mutual Assessment Treatment
sum  Formidability_z maxgrip maxchest flexedbicep perceivedstrength age height weight bmi risk if treatment==1& period==uniqueperiod

*Balance tests and CIs
preserve
collapse Formidability_z maxgrip maxchest flexedbicep perceivedstrength age height weight bmi  risk treatment, by(participantid)
ttest age==0
ttest Formidability_z==0
ttest height==0
ttest weight==0
ttest bmi==0
ttest maxgrip==0
ttest maxchest==0
ttest flexedbicep==0
ttest perceivedstrength==0
ttest risk==0
restore 

**********************************************************
***********Table S3: Negative binomial regressions for the duration of contests (Mutual Assessment Treatment in Study 3) ***********
**********************************************************
which regsave //user written stata command for saving regression output
*Prepare a matrix for 99-percent CI*
set more off
matrix alphaCI = J(4,8,.)
matrix colnames alphaCI =  model1 model2 model3 model4 model5 model6 model7 model8
matrix rownames alphaCI = alpha alphase lower upper

preserve 
keep if treatment==1 & contestantincluded==1
set level 99
*model1
vce2way nbreg finalcosts absdstrength per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 per11 per12 per13  , cluster(id1 id2)
eststo model1
matrix tempmat=r(table)
local colno=colsof(tempmat)
matrix alphaCI[1,1]=tempmat[1,`colno']
matrix alphaCI[2,1]=tempmat[2,`colno']
matrix alphaCI[3,1]=tempmat[5,`colno']
matrix alphaCI[4,1]=tempmat[6,`colno']
*model2
vce2way nbreg finalcosts absdheight per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 per11 per12 per13  , cluster(id1 id2)
eststo model2
matrix tempmat=r(table)
local colno=colsof(tempmat)
matrix alphaCI[1,2]=tempmat[1,`colno']
matrix alphaCI[2,2]=tempmat[2,`colno']
matrix alphaCI[3,2]=tempmat[5,`colno']
matrix alphaCI[4,2]=tempmat[6,`colno']
*model3
vce2way nbreg finalcosts absdstrength absdheight c.absdstrength#c.absdheight per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 per11 per12 per13  , cluster(id1 id2)
eststo model3
matrix tempmat=r(table)
local colno=colsof(tempmat)
matrix alphaCI[1,3]=tempmat[1,`colno']
matrix alphaCI[2,3]=tempmat[2,`colno']
matrix alphaCI[3,3]=tempmat[5,`colno']
matrix alphaCI[4,3]=tempmat[6,`colno']
*model4
vce2way nbreg finalcosts absdbmi per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 per11 per12 per13  , cluster(id1 id2)
eststo model4
matrix tempmat=r(table)
local colno=colsof(tempmat)
matrix alphaCI[1,4]=tempmat[1,`colno']
matrix alphaCI[2,4]=tempmat[2,`colno']
matrix alphaCI[3,4]=tempmat[5,`colno']
matrix alphaCI[4,4]=tempmat[6,`colno']
*model5
vce2way nbreg finalcosts absdrisk per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 per11 per12 per13  , cluster(id1 id2)
eststo model5
matrix tempmat=r(table)
local colno=colsof(tempmat)
matrix alphaCI[1,5]=tempmat[1,`colno']
matrix alphaCI[2,5]=tempmat[2,`colno']
matrix alphaCI[3,5]=tempmat[5,`colno']
matrix alphaCI[4,5]=tempmat[6,`colno']
*model6
vce2way nbreg finalcosts absdstrength absdrisk per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 per11 per12 per13  , cluster(id1 id2)
eststo model6
matrix tempmat=r(table)
local colno=colsof(tempmat)
matrix alphaCI[1,6]=tempmat[1,`colno']
matrix alphaCI[2,6]=tempmat[2,`colno']
matrix alphaCI[3,6]=tempmat[5,`colno']
matrix alphaCI[4,6]=tempmat[6,`colno']
**model7: leave out as minor detail beyond models 5 and 8
*vce2way nbreg finalcosts strongrisk weakrisk per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 per11 per12 per13  , cluster(id1 id2)
*eststo model7
*matrix tempmat=r(table)
*local colno=colsof(tempmat)
*matrix alphaCI[1,7]=tempmat[1,`colno']
*matrix alphaCI[2,7]=tempmat[2,`colno']
*matrix alphaCI[3,7]=tempmat[5,`colno']
*matrix alphaCI[4,7]=tempmat[6,`colno']
*model7
vce2way nbreg finalcosts strongformid weakformid per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 per11 per12 per13  , cluster(id1 id2)
eststo model7
matrix tempmat=r(table)
local colno=colsof(tempmat)
matrix alphaCI[1,7]=tempmat[1,`colno']
matrix alphaCI[2,7]=tempmat[2,`colno']
matrix alphaCI[3,7]=tempmat[5,`colno']
matrix alphaCI[4,7]=tempmat[6,`colno']
*model8
vce2way nbreg finalcosts strongrisk weakrisk strongformid weakformid taller_z shorter_z per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 per11 per12 per13  , cluster(id1 id2)
eststo model8
matrix tempmat=r(table)
local colno=colsof(tempmat)
matrix alphaCI[1,8]=tempmat[1,`colno']
matrix alphaCI[2,8]=tempmat[2,`colno']
matrix alphaCI[3,8]=tempmat[5,`colno']
matrix alphaCI[4,8]=tempmat[6,`colno']

restore
matrix list alphaCI
esttab using duration_MAT.rtf, replace /// 
	addnote("Robust standard errors") drop(per*) label nonumber mtitle("Model 1" "Model 2" "Model 3" "Model 4" "Model 5" "Model 6" "Model 7" "Model 8") b(3) se(3) se starlevels(* 0.1 ** 0.05 *** 0.01)  
*note: 95 clusters in id 1 and 96 clusters in id2 (one person had to be dropped, hence asymmetric) because there are 7 sessions and the max and min subject ids from each session only appear once, either in id1 or id2, whereas all other subject ids appear in both id1 and id2
eststo clear

**********************************************************
*model selection
**********************************************************
*Note: Can't use cluster with LR test as the likelihood does not account for the clustering, alternative is to use Wald test, 
*ie. CI reported for alpha, level set to 1% (see above)

*LR test
*model1
poisson finalcosts absdstrength per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 per11 per12 per13 if treatment==1 & contestantincluded==1 , robust
est store poisson
nbreg finalcosts absdstrength per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 per11 per12 per13 if treatment==1 & contestantincluded==1 , robust
est store nbreg
lrtest poisson nbreg, stats force
eststo clear

*model2
poisson finalcosts absdheight per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 per11 per12 per13  if treatment==1 & contestantincluded==1, robust
est store poisson
nbreg finalcosts absdheight per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 per11 per12 per13  if treatment==1 & contestantincluded==1, robust
est store nbreg
lrtest poisson nbreg, stats force
eststo clear

*model3
poisson finalcosts absdstrength absdheight c.absdstrength#c.absdheight per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 per11 per12 per13  if treatment==1 & contestantincluded==1, robust
est store poisson
nbreg finalcosts absdstrength absdheight c.absdstrength#c.absdheight per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 per11 per12 per13  if treatment==1 & contestantincluded==1, robust
est store nbreg
lrtest poisson nbreg, stats force
eststo clear

*model4
poisson finalcosts absdbmi per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 per11 per12 per13 if treatment==1 & contestantincluded==1 , robust
est store poisson
nbreg finalcosts absdbmi per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 per11 per12 per13 if treatment==1 & contestantincluded==1 , robust
est store nbreg
lrtest poisson nbreg, stats force
eststo clear


*model5
poisson finalcosts absdrisk per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 per11 per12 per13  if treatment==1 & contestantincluded==1, robust
est store poisson
nbreg finalcosts absdrisk per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 per11 per12 per13  if treatment==1 & contestantincluded==1, robust
est store nbreg
lrtest poisson nbreg, stats force
eststo clear

*model6
poisson finalcosts absdstrength absdrisk per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 per11 per12 per13  if treatment==1 & contestantincluded==1, robust
est store poisson
nbreg finalcosts absdstrength absdrisk per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 per11 per12 per13  if treatment==1 & contestantincluded==1, robust
est store nbreg
lrtest poisson nbreg, stats force
eststo clear

*model7
poisson finalcosts strongformid weakformid per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 per11 per12 per13  if treatment==1 & contestantincluded==1, robust
est store poisson
nbreg finalcosts strongformid weakformid per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 per11 per12 per13  if treatment==1 & contestantincluded==1, robust
est store nbreg
lrtest poisson nbreg, stats force
eststo clear

*model8
poisson finalcosts strongrisk weakrisk strongformid weakformid taller_z shorter_z per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 per11 per12 per13  if treatment==1 & contestantincluded==1, robust
est store poisson
nbreg finalcosts strongrisk weakrisk strongformid weakformid taller_z shorter_z per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 per11 per12 per13  if treatment==1 & contestantincluded==1, robust
est store nbreg
lrtest poisson nbreg, stats force
eststo clear

*CONCLUSION: alpha significantly differs from 0 -> overdispersion, reject poisson in favor of neg.binomial model


************************************************************************
***********Table S5: Logistic regressions on whether the strongest contestant wins (Mutual Assessment Treatment in Study 3)**********
************************************************************************


preserve
keep if treatment==1 & contestantincluded==1 & draw!=1 
vce2way logit strongwins absdstrength per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 per11 per12 per13,  cluster(id1 id2)
eststo model1
vce2way logit strongwins absdstrength absdstrength2 per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 per11 per12 per13, cluster(id1 id2)
eststo model2
//To test that absdstrength and absdstrength are jointly 0:
test absdstrength absdstrength2

vce2way logit strongwins strongrisk weakrisk per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 per11 per12 per13, cluster(id1 id2)
eststo model3
vce2way logit strongwins absdstrength absdstrength2 strongrisk weakrisk per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 per11 per12 per13, cluster(id1 id2)
eststo model4
keep if treatment==1 & contestantincluded==1 & finalcosts==0 
vce2way logit strongwins absdstrength per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 per11 per12 per13,  cluster(id1 id2)
eststo model5
vce2way logit strongwins strongrisk weakrisk per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 per11 per12 per13, cluster(id1 id2)
eststo model6
vce2way logit strongwins absdstrength strongrisk weakrisk per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 per11 per12 per13, cluster(id1 id2)
eststo model7
restore

esttab using winner_strength_MAT.rtf, replace ///
	addnote("Logistic regression") drop(per*) label nonumber mtitle("Model 1" "Model 2" "Model 3" "Model 4") b(3) se(3) se starlevels(* 0.10 ** 0.05 *** 0.01)
eststo clear





***********************************************************************
***********Table S7: (Study 3) Logistic regressions on whether the tallest contestant wins (Mutual Assessment Treatment)****************
***********************************************************************

*exclude cases where sameheight==1
by gameid, sort: egen sumtallest = sum(tallest) 

preserve
keep if tallest==1 & sumtallest==1 & treatment==1 & draw!=1

vce2way logit tallwin absdheight per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 per11 per12 per13 , cluster(id1 id2)
eststo model1
vce2way logit tallwin absdheight absdheight2 per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 per11 per12 per13 , cluster(id1 id2)
eststo model2
vce2way logit tallwin tallrisk shortrisk per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 per11 per12 per13 , cluster(id1 id2)
eststo model3
vce2way logit tallwin absdheight tallrisk shortrisk per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 per11 per12 per13 , cluster(id1 id2)
eststo model4
restore 
esttab using winner_height_MAT.rtf, replace ///
	addnote("Logistic regression") drop(per*) label nonumber mtitle("Model 1" "Model 2" "Model 3" "Model 4") b(3) se(3) se starlevels(* 0.10 ** 0.05 *** 0.01) 
eststo clear

************************************************************************
***********Table S8: (Study 3) Logistic regressions on whether the leanest contestant wins (Mutual Assessment Treatment)********************
************************************************************************

*exclude cases where samebmi==1
by gameid, sort: egen sumleanest = sum(leanest)

preserve
keep if treatment==1 & leanest==1 & sumleanest==1 & draw!=1

vce2way logit leanwin absdbmi per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 per11 per12 per13 , cluster(id1 id2)
eststo model1
vce2way logit leanwin absdbmi absdbmi2 per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 per11 per12 per13 , cluster(id1 id2)
eststo model2
vce2way logit leanwin leanrisk heavyrisk per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 per11 per12 per13 , cluster(id1 id2)
eststo model3
vce2way logit leanwin absdbmi leanrisk heavyrisk per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 per11 per12 per13 , cluster(id1 id2)
eststo model4
restore 
esttab using winner_bmi_MAT.rtf,  replace ///
	addnote("Logistic regression") drop(per*) label nonumber mtitle("Model 1" "Model 2" "Model 3" "Model 4") b(3) se(3) se starlevels(* 0.10 ** 0.05 *** 0.01) 
eststo clear









******************************************************************************************************************************************************************************
******************************************************************************************************************************************************************************

**********************************************************
***********Study 4***********
**********Note: results from the replication study (Study 3) are stated below, after the results from Study 4 **********
**********************************************************
cd "..\"

use data_study4,clear
set more off

cd "./results"

*Note: in the following double precision is important so that c== actually works, because RHS varable is double precision, whereas standard generate creates single precision and rounding causes == to fail
set type double


***************************************************************
***********Figure 5: (Study 4) Mediation***********
***************************************************************
which combomarginsplot
preserve
keep if contestantincluded==1
vce2way regress absdaggression absdstrength per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 , cluster(id1 id2)
predict double stragg_res, residuals

quietly: vce2way nbreg finalcosts absdstrength absdaggression per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 , cluster(id1 id2)
quietly: margins, at(absdstrength=(0(0.1)3.9)) atmeans saving(fileM1, replace)

quietly: vce2way nbreg finalcosts absdstrength stragg_res per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 , cluster(id1 id2)
quietly: margins, at(absdstrength=(0(0.1)3.9)) atmeans saving(fileM2, replace)
combomarginsplot fileM2 fileM1 , labels("Total effect" "Direct effect")  file1opts(recast(line) lw(vthick) lcolor(gray) lpattern(dash)) file2opts(recast(line) lw(vthick) lcolor(black)) noci legend(colfirst) legend(size(medium) ring(0) position(2)) ytitle("Predicted duration in seconds") xtitle("Absolute difference in strength") title("") xscale(titlegap(3)) yscale(titlegap(2)) scheme(lean1) saving(filefig4_3, replace)
restore




********************************************************************************************************************
***** Table S9: Descriptive statistics for the personality feature ratings in Study 4 (unstandardized)*****
********************************************************************************************************************
sum mean_aggression mean_getway mean_friend mean_intelligent if obsno==1

********************************************************************************************************************
***** Table S10: (Study 4) Correlations of average personality feature ratings among each other and with strength (standardized) *****
********************************************************************************************************************
pwcorr  z_mean_aggression z_mean_getway z_mean_friend z_mean_intelligent Formidability_z if obsno==1 , sig obs




********************************************************************************************************************
*Table S9 Potential mediators of the effect of differences in strength on contest duration in the Mutual Assessment Treatment (Studies 2 & 4)*
********************************************************************************************************************
which matsave
* if not installed, use findit matsave and load matsave from http://fmwww.bc.edu/RePEc/bocode/m
** WARNING: If the path "dir\matrices\" does not exist, matsave will result in the error **

*Note: this report results based on the pre-saved bootstrap samples
*To run the bootstrap use the code below (search for "RE-RUN Bootstrap")*

cd "..\"


*Prepare a matrix with coefficients, bootstrap standard errors and bias corrected confidence intervals (levels 90,95,99)*
matrix mediation = J(48,6,.)
*								1	2		3		4		5				6						
matrix colnames mediation =  total direct indirect perc indirectisentangle mediationperc  
*								1		2	3	4	5	 6		7	8		
matrix rownames mediation = aggression se ll90 ul90 ll95 ul95 ll99 ul99 ///
getway se ll90 ul90 ll95 ul95 ll99 ul99 friend se ll90 ul90 ll95 ul95 ll99 ul99 ///
intelligent  se ll90 ul90 ll95 ul95 ll99 ul99 ///
totalindir se ll90 ul90 ll95 ul95 ll99 ul99 direct se ll90 ul90 ll95 ul95 ll99 ul99


*** absdaggression*****
quiet{
	*The estimate in our sample
	preserve
	keep if contestantincluded==1
	
	quietly: vce2way regress absdaggression absdstrength per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 , cluster(id1 id2)
	quietly: predict double str_res, residuals

	vce2way nbreg finalcosts absdstrength absdaggression per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 , cluster(id1 id2)
	matrix erase = e(b)
	matrix khbstats = J(1,3,.)
	matrix khbstats[1,2]=erase[1,1]

	vce2way nbreg finalcosts absdstrength str_res per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 , cluster(id1 id2)
	matrix erase = e(b)
	matrix khbstats[1,1]=erase[1,1]
	matrix khbstats[1,3]=khbstats[1,1]-khbstats[1,2]
	matrix drop erase
	matrix list khbstats
	matrix mediation[1,1]=khbstats
	* WARNING: If the path "dir\matrices\" does not exist, matsave will result in the error
	set more off
	matsave mediation, replace p("matrices") dropall
	restore

	set more off
	matload mediation,  p("matrices")  saving over
	
	*Bootstrap se and bias corrected CIs (recognizes what to compute from "estimate in our sample" template above)
	bstat using khbdagg, stat(khbstats) n(205)  level(90)
	estat bootstrap, all 
	
	matrix mediation[2,1]=e(se)
	matrix mediation[3,1]=e(ci_bc)
		
	bstat using khbdagg, stat(khbstats) n(205)  level(95)
	estat bootstrap, all 
	matrix mediation[5,1]=e(ci_bc)
	
	bstat using khbdagg, stat(khbstats) n(205)  level(99)
	estat bootstrap, all 
	matrix mediation[7,1]=e(ci_bc)
}

*** absdgetway*****	
quiet{
local counter=8
	
	*The estimate in our sample
	preserve
	keep if contestantincluded==1
	
	quietly: vce2way regress absdgetway absdstrength per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 , cluster(id1 id2)
	quietly: predict double str_res, residuals

	vce2way nbreg finalcosts absdstrength absdgetway per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 , cluster(id1 id2)
	matrix erase = e(b)
	matrix khbstats = J(1,3,.)
	matrix khbstats[1,2]=erase[1,1]

	vce2way nbreg finalcosts absdstrength str_res per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 , cluster(id1 id2)
	matrix erase = e(b)
	matrix khbstats[1,1]=erase[1,1]
	matrix khbstats[1,3]=khbstats[1,1]-khbstats[1,2]
	matrix drop erase
	matrix list khbstats
	matrix mediation[`counter'+1,1]=khbstats
	* WARNING: If the path "dir\matrices\" does not exist, matsave will result in the error
	set more off
	matsave mediation, replace p("matrices") dropall
	restore

	set more off
	matload mediation,  p("matrices")  saving over
	
	*Bootstrap se and bias corrected CIs (recognizes what to compute from "estimate in our sample" template above)
	bstat using khbdget, stat(khbstats) n(205)  level(90)
	estat bootstrap, all 
	
	matrix mediation[`counter'+2,1]=e(se)
	matrix mediation[`counter'+3,1]=e(ci_bc)
		
	bstat using khbdget, stat(khbstats) n(205)  level(95)
	estat bootstrap, all 
	matrix mediation[`counter'+5,1]=e(ci_bc)
	
	bstat using khbdget, stat(khbstats) n(205)  level(99)
	estat bootstrap, all 
	matrix mediation[`counter'+7,1]=e(ci_bc)	
}


*** absdfriend*****	
quiet{
local counter=16
	
	*The estimate in our sample
	preserve
	keep if contestantincluded==1
	
	quietly: vce2way regress absdfriend absdstrength per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 , cluster(id1 id2)
	quietly: predict double str_res, residuals

	vce2way nbreg finalcosts absdstrength absdfriend per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 , cluster(id1 id2)
	matrix erase = e(b)
	matrix khbstats = J(1,3,.)
	matrix khbstats[1,2]=erase[1,1]

	vce2way nbreg finalcosts absdstrength str_res per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 , cluster(id1 id2)
	matrix erase = e(b)
	matrix khbstats[1,1]=erase[1,1]
	matrix khbstats[1,3]=khbstats[1,1]-khbstats[1,2]
	matrix drop erase
	matrix list khbstats
	matrix mediation[`counter'+1,1]=khbstats
	* WARNING: If the path "dir\matrices\" does not exist, matsave will result in the error
	set more off
	matsave mediation, replace p("matrices") dropall
	restore

	set more off
	matload mediation,  p("matrices")  saving over
	
	*Bootstrap se and bias corrected CIs (recognizes what to compute from "estimate in our sample" template above)
	bstat using khbdfri, stat(khbstats) n(205)  level(90)
	estat bootstrap, all 
	
	matrix mediation[`counter'+2,1]=e(se)
	matrix mediation[`counter'+3,1]=e(ci_bc)
		
	bstat using khbdfri, stat(khbstats) n(205)  level(95)
	estat bootstrap, all 
	matrix mediation[`counter'+5,1]=e(ci_bc)
	
	bstat using khbdfri, stat(khbstats) n(205)  level(99)
	estat bootstrap, all 
	matrix mediation[`counter'+7,1]=e(ci_bc)
}	

*** absdintelligent*****	
quiet{
local counter=24
	
	*The estimate in our sample
	preserve
	keep if contestantincluded==1
	
	quietly: vce2way regress absdintelligent absdstrength per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 , cluster(id1 id2)
	quietly: predict double str_res, residuals

	vce2way nbreg finalcosts absdstrength absdintelligent per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 , cluster(id1 id2)
	matrix erase = e(b)
	matrix khbstats = J(1,3,.)
	matrix khbstats[1,2]=erase[1,1]

	vce2way nbreg finalcosts absdstrength str_res per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 , cluster(id1 id2)
	matrix erase = e(b)
	matrix khbstats[1,1]=erase[1,1]
	matrix khbstats[1,3]=khbstats[1,1]-khbstats[1,2]
	matrix drop erase
	matrix list khbstats
	matrix mediation[`counter'+1,1]=khbstats
	* WARNING: If the path "dir\matrices\" does not exist, matsave will result in the error
	set more off
	matsave mediation, replace p("matrices") dropall
	restore

	set more off
	matload mediation,  p("matrices")  saving over
	
	*Bootstrap se and bias corrected CIs (recognizes what to compute from "estimate in our sample" template above)
	bstat using khbdint, stat(khbstats) n(205)  level(90)
	estat bootstrap, all 
	
	matrix mediation[`counter'+2,1]=e(se)
	matrix mediation[`counter'+3,1]=e(ci_bc)
		
	bstat using khbdint, stat(khbstats) n(205)  level(95)
	estat bootstrap, all 
	matrix mediation[`counter'+5,1]=e(ci_bc)
	
	bstat using khbdint, stat(khbstats) n(205)  level(99)
	estat bootstrap, all 
	matrix mediation[`counter'+7,1]=e(ci_bc)

*** mediation precentage ***
local j=1
forvalues i=1(1)4{
	matrix mediation[`j',4]=mediation[`j',3]/mediation[`j',1]*100
	local `j'=`j'+8
	}

	matrix mediation[1,4]=mediation[1,3]/mediation[1,1]*100
	matrix mediation[9,4]=mediation[9,3]/mediation[9,1]*100
	matrix mediation[17,4]=mediation[17,3]/mediation[17,1]*100
	matrix mediation[25,4]=mediation[25,3]/mediation[25,1]*100
}


*** disentangled *****	
quiet{
	*The estimate in our sample
	preserve
	keep if contestantincluded==1
	
	matrix khbstats = J(5,7,.)
	matrix colnames khbstats =  zonx  coeff indirect direct total percindirect perctotal
	matrix rownames khbstats = aggression getway friend intelligent sum


	*effect of x on mediator z
	regress absdaggression absdstrength per1 per2 per3 per4 per5 per6 per7 per8 per9 per10
	matrix erase = e(b)
	matrix khbstats[1,1]=erase[1,1]
	matrix drop erase
	quietly: predict double stragg_res, residuals

	regress absdgetway absdstrength  per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 
	matrix erase = e(b)
	matrix khbstats[2,1]=erase[1,1]
	matrix drop erase
	quietly: predict double strget_res, residuals

	regress absdfriend absdstrength  per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 
	matrix erase = e(b)
	matrix khbstats[3,1]=erase[1,1]
	matrix drop erase
	quietly: predict double strfriend_res, residuals

	regress  absdintelligent absdstrength per1 per2 per3 per4 per5 per6 per7 per8 per9 per10
	matrix erase = e(b)
	matrix khbstats[4,1]=erase[1,1]
	matrix drop erase
	quietly: predict double strintell_res, residuals

	*coefficient for computation of indirect effect
	vce2way nbreg finalcosts absdstrength absdaggression absdgetway absdfriend absdintelligent  per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 , cluster(id1 id2)
	matrix erase = e(b)
	matrix khbstats[1,2]=erase[1,2]
	matrix khbstats[2,2]=erase[1,3]
	matrix khbstats[3,2]=erase[1,4]
	matrix khbstats[4,2]=erase[1,5]

	*indirect effect
	matrix khbstats[1,3]=erase[1,2]*khbstats[1,1]
	matrix khbstats[2,3]=erase[1,3]*khbstats[2,1]
	matrix khbstats[3,3]=erase[1,4]*khbstats[3,1]
	matrix khbstats[4,3]=erase[1,5]*khbstats[4,1]
	matrix khbstats[5,3]=khbstats[1,3]+khbstats[2,3]+khbstats[3,3]+khbstats[4,3]

	*direct effect
	forvalues i=1(1)4 {
	matrix khbstats[`i',4]=erase[1,1]
	}

	*total effect=direct effect + sum of indirect effects
	forvalues i=1(1)4 {
	matrix khbstats[`i',5]=khbstats[1,4]+khbstats[5,3]
	}

	*the contribution percentage of each mediator to the indirect eﬀect, sums up to 100
	forvalues i=1(1)4 {
	matrix khbstats[`i',6]=khbstats[`i',3]/khbstats[5,3]*100
	}
	matrix khbstats[5,6]=khbstats[1,6]+khbstats[2,6]+khbstats[3,6]+khbstats[4,6]

	*the mediation percentage (contribution percentage of each mediator to the total eﬀect, sums up to the overall confounding percentage=(sum of indirect effects)/(total effect))
	matrix khbstats[5,7]=khbstats[5,3]/khbstats[1,5]*100 //overall confounding percentage
	forvalues i=1(1)4 {
	matrix khbstats[`i',7]=khbstats[`i',6]*khbstats[5,7]/100
	}
	matrix drop erase
	matrix list khbstats
	*record results to return 
	*order: (1-4:indir effects, 7-9: mediation percentage) 1.absdaggression, 2.absdgetway, 3.absdfriend, 4.absdintelligence 5.total indirect effects 6.direct effect of absdstrength 7.absdaggression, 8.absdgetway, 9.absdfriend, 10.absdintelligence 11.total
	matrix disentangled = [khbstats[1,3],khbstats[2,3],khbstats[3,3],khbstats[4,3],khbstats[5,3],khbstats[1,4],khbstats[1,7],khbstats[2,7],khbstats[3,7],khbstats[4,7],khbstats[5,7]]
	matrix list disentangled
	

	matrix mediation[1,5]=disentangled[1,1]
	matrix mediation[9,5]=disentangled[1,2]
	matrix mediation[17,5]=disentangled[1,3]
	matrix mediation[25,5]=disentangled[1,4]
	matrix mediation[33,5]=disentangled[1,5]
	matrix mediation[41,5]=disentangled[1,6]
	
	matrix mediation[1,6]=disentangled[1,7]
	matrix mediation[9,6]=disentangled[1,8]
	matrix mediation[17,6]=disentangled[1,9]
	matrix mediation[25,6]=disentangled[1,10]
	matrix mediation[33,6]=disentangled[1,11]
		
	* WARNING: If the path "dir\matrices\" does not exist, matsave will result in the error
	set more off
	matsave mediation, replace p("matrices") dropall
	* WARNING: If the path "dir\matrices\" does not exist, matsave will result in the error
	set more off
	matsave mediation, replace p("matrices") dropall
	restore

	set more off
	matload mediation,  p("matrices")  saving over
	
	*Bootstrap se and bias corrected CIs (recognizes what to compute from "estimate in our sample" template above)
	bstat using khbdisentangle, stat(disentangled) n(205) level(90)
	estat bootstrap, all
	matrix tempse=e(se)
	matrix tempci=e(ci_bc)

	local j=2
	forvalues i=1(1)6{
		local j1=`j'+1
		local j2=`j'+2	
		matrix mediation[`j',5]=tempse[1,`i']
		matrix mediation[`j1',5]=tempci[1,`i']
		matrix mediation[`j2',5]=tempci[2,`i']
		local j=`j'+8
		}
	bstat using khbdisentangle, stat(disentangled) n(205) level(95)
	estat bootstrap, all
	matrix tempci=e(ci_bc)

	local j=5
	forvalues i=1(1)6{
		local j1=`j'+1
		matrix mediation[`j',5]=tempci[1,`i']
		matrix mediation[`j1',5]=tempci[2,`i']
		local j=`j'+8
		}
		
	bstat using khbdisentangle, stat(disentangled) n(205) level(99)
	estat bootstrap, all
	matrix tempci=e(ci_bc)

	local j=7
	forvalues i=1(1)6{
		local j1=`j'+1
		matrix mediation[`j',5]=tempci[1,`i']
		matrix mediation[`j1',5]=tempci[2,`i']
		local j=`j'+8
		}
}
*note: missing values in matrix are -999		
matrix list mediation
		
		
*****************************************************************************************************************************
** RE-RUN Bootstrap **
** Run separately to generate bootstrap sample data for table **
*****************************************************************************************************************************


** Write bootstrap samples into results folder (ie don't overwrite original data)**
cd "./results"

*** absdaggression*****
quiet{
capture program drop myboot2
program define myboot2, rclass
	preserve 
		tsset participantid period
		*sample from participant panel observations to "fill" sessions with resampling
		bsample, cluster(participantid) idcluster(newparticid) strata(session)
		*assign new gameid for duplicate games
		by gameid, sort: gen seqrepetition = _n
		egen newgameid=group(seqrepetition gameid)
		*(Not needed for bootstrap) store in each row both subject ids for each game, for clustering
		by newgameid, sort: egen newid1=min(newparticid)
		by newgameid, sort: egen newid2=max(newparticid)
		*generate dummy for which contests to include: drop 205 rows to obtain 205 contests, balanced across periods
		gen randno=runiform()
		sort period randno
		by period: gen pno=_n
		by period: gen pmax=_N/2
		drop if pno>pmax
		drop pno pmax
 
		 regress absdaggression absdstrength per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 
		predict double bootstragg_res, residuals

		nbreg finalcosts absdstrength absdaggression per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 
		matrix erase = e(b)
		matrix khbstats = J(1,3,.)
		matrix khbstats[1,2]=erase[1,1]

		nbreg finalcosts absdstrength bootstragg_res per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 
		matrix erase = e(b)
		matrix khbstats[1,1]=erase[1,1]
		matrix khbstats[1,3]=khbstats[1,1]-khbstats[1,2]
		matrix drop erase
		matrix erase = e(b)
		return scalar khbred = khbstats[1,1]
		return scalar khbfull = khbstats[1,2]
		return scalar khbdif = khbstats[1,3]
	restore
end
}
*Run bootstrap
preserve
	simulate khbred=r(khbred) khbfull=r(khbfull) khbdif=r(khbdif), ///
	reps(1000) seed(12345) saving(khbdagg, replace): myboot2
restore

*** absdgetway*****
quiet{
capture program drop myboot2
program define myboot2, rclass
	preserve 
		tsset participantid period
		*sample from participant panel observations to "fill" sessions with resampling
		bsample, cluster(participantid) idcluster(newparticid) strata(session)
		*assign new gameid for duplicate games
		by gameid, sort: gen seqrepetition = _n
		egen newgameid=group(seqrepetition gameid)
		*(Not needed for bootstrap) store in each row both subject ids for each game, for clustering
		by newgameid, sort: egen newid1=min(newparticid)
		by newgameid, sort: egen newid2=max(newparticid)
		*generate dummy for which contests to include: drop 205 rows to obtain 205 contests, balanced across periods
		gen randno=runiform()
		sort period randno
		by period: gen pno=_n
		by period: gen pmax=_N/2
		drop if pno>pmax
		drop pno pmax
		  
		 regress absdgetway absdstrength per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 
		predict double bootstragg_res, residuals

		nbreg finalcosts absdstrength absdgetway per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 
		matrix erase = e(b)
		matrix khbstats = J(1,3,.)
		matrix khbstats[1,2]=erase[1,1]

		nbreg finalcosts absdstrength bootstragg_res per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 
		matrix erase = e(b)
		matrix khbstats[1,1]=erase[1,1]
		matrix khbstats[1,3]=khbstats[1,1]-khbstats[1,2]
		matrix drop erase
		matrix erase = e(b)
		return scalar khbred = khbstats[1,1]
		return scalar khbfull = khbstats[1,2]
		return scalar khbdif = khbstats[1,3]
	restore
end
}
*Run bootstrap
preserve
	simulate khbred=r(khbred) khbfull=r(khbfull) khbdif=r(khbdif), ///
	reps(1000) seed(12345) saving(khbdget, replace): myboot2
restore
      
*** absdfriend*****
quiet{
capture program drop myboot2
program define myboot2, rclass
	preserve 
		tsset participantid period
		*sample from participant panel observations to "fill" sessions with resampling
		bsample, cluster(participantid) idcluster(newparticid) strata(session)
		*assign new gameid for duplicate games
		by gameid, sort: gen seqrepetition = _n
		egen newgameid=group(seqrepetition gameid)
		*(Not needed for bootstrap) store in each row both subject ids for each game, for clustering
		by newgameid, sort: egen newid1=min(newparticid)
		by newgameid, sort: egen newid2=max(newparticid)
		*generate dummy for which contests to include: drop 205 rows to obtain 205 contests, balanced across periods
		gen randno=runiform()
		sort period randno
		by period: gen pno=_n
		by period: gen pmax=_N/2
		drop if pno>pmax
		drop pno pmax
		  
		 regress absdfriend absdstrength per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 
		predict double bootstragg_res, residuals

		nbreg finalcosts absdstrength absdfriend per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 
		matrix erase = e(b)
		matrix khbstats = J(1,3,.)
		matrix khbstats[1,2]=erase[1,1]

		nbreg finalcosts absdstrength bootstragg_res per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 
		matrix erase = e(b)
		matrix khbstats[1,1]=erase[1,1]
		matrix khbstats[1,3]=khbstats[1,1]-khbstats[1,2]
		matrix drop erase
		matrix erase = e(b)
		return scalar khbred = khbstats[1,1]
		return scalar khbfull = khbstats[1,2]
		return scalar khbdif = khbstats[1,3]
	restore
end
}
*Run bootstrap
preserve
	simulate khbred=r(khbred) khbfull=r(khbfull) khbdif=r(khbdif), ///
	reps(1000) seed(12345) saving(khbdfri, replace): myboot2
restore

*** absdintelligent*****
quiet{
capture program drop myboot2
program define myboot2, rclass
	preserve 
		tsset participantid period
		*sample from participant panel observations to "fill" sessions with resampling
		bsample, cluster(participantid) idcluster(newparticid) strata(session)
		*assign new gameid for duplicate games
		by gameid, sort: gen seqrepetition = _n
		egen newgameid=group(seqrepetition gameid)
		*(Not needed for bootstrap) store in each row both subject ids for each game, for clustering
		by newgameid, sort: egen newid1=min(newparticid)
		by newgameid, sort: egen newid2=max(newparticid)
		*generate dummy for which contests to include: drop 205 rows to obtain 205 contests, balanced across periods
		gen randno=runiform()
		sort period randno
		by period: gen pno=_n
		by period: gen pmax=_N/2
		drop if pno>pmax
		drop pno pmax
		  
		 regress absdintelligent absdstrength per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 
		predict double bootstragg_res, residuals

		nbreg finalcosts absdstrength absdintelligent per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 
		matrix erase = e(b)
		matrix khbstats = J(1,3,.)
		matrix khbstats[1,2]=erase[1,1]

		nbreg finalcosts absdstrength bootstragg_res per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 
		matrix erase = e(b)
		matrix khbstats[1,1]=erase[1,1]
		matrix khbstats[1,3]=khbstats[1,1]-khbstats[1,2]
		matrix drop erase
		matrix erase = e(b)
		return scalar khbred = khbstats[1,1]
		return scalar khbfull = khbstats[1,2]
		return scalar khbdif = khbstats[1,3]
	restore
end
}
*Run bootstrap
preserve
	simulate khbred=r(khbred) khbfull=r(khbfull) khbdif=r(khbdif), ///
	reps(1000) seed(12345) saving(khbdint, replace): myboot2
restore

*** all mediators at once, disentangle contributions to mediation*****
use data_study3,clear
set more off


quiet{
capture program drop myboot2
program define myboot2, rclass
	preserve 

		matrix khbstats = J(5,7,.)
		matrix colnames khbstats =  zonx  coeff indirect direct total percindirect perctotal
		matrix rownames khbstats = aggression getway friend intelligent sum

		tsset participantid period
		*sample from participant panel observations to "fill" sessions with resampling
		bsample, cluster(participantid) idcluster(newparticid) strata(session)
		*assign new gameid for duplicate games
		by gameid, sort: gen seqrepetition = _n
		egen newgameid=group(seqrepetition gameid)
		*(Not needed for bootstrap) store in each row both subject ids for each game, for clustering
		by newgameid, sort: egen newid1=min(newparticid)
		by newgameid, sort: egen newid2=max(newparticid)
		*generate dummy for which contests to include: drop 205 rows to obtain 205 contests, balanced across periods
		gen randno=runiform()
		sort period randno
		by period: gen pno=_n
		by period: gen pmax=_N/2
		drop if pno>pmax
		drop pno pmax


		*effect of x on mediator z
		regress absdaggression absdstrength per1 per2 per3 per4 per5 per6 per7 per8 per9 per10
		matrix erase = e(b)
		matrix khbstats[1,1]=erase[1,1]
		matrix drop erase
		predict double stragg_res, residuals

		regress absdgetway absdstrength  per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 
		matrix erase = e(b)
		matrix khbstats[2,1]=erase[1,1]
		matrix drop erase
		predict double strget_res, residuals

		regress absdfriend absdstrength  per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 
		matrix erase = e(b)
		matrix khbstats[3,1]=erase[1,1]
		matrix drop erase
		predict double strfriend_res, residuals

		regress  absdintelligent absdstrength per1 per2 per3 per4 per5 per6 per7 per8 per9 per10
		matrix erase = e(b)
		matrix khbstats[4,1]=erase[1,1]
		matrix drop erase
		predict double strintell_res, residuals

		*coefficient for computation of indirect effect
		nbreg finalcosts absdstrength absdaggression absdgetway absdfriend absdintelligent  per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 
		matrix erase = e(b)
		matrix khbstats[1,2]=erase[1,2]
		matrix khbstats[2,2]=erase[1,3]
		matrix khbstats[3,2]=erase[1,4]
		matrix khbstats[4,2]=erase[1,5]


		*indirect effect
		matrix khbstats[1,3]=erase[1,2]*khbstats[1,1]
		matrix khbstats[2,3]=erase[1,3]*khbstats[2,1]
		matrix khbstats[3,3]=erase[1,4]*khbstats[3,1]
		matrix khbstats[4,3]=erase[1,5]*khbstats[4,1]
		matrix khbstats[5,3]=khbstats[1,3]+khbstats[2,3]+khbstats[3,3]+khbstats[4,3]

		*direct effect
		forvalues i=1(1)4 {
		matrix khbstats[`i',4]=erase[1,1]
		}

		*total effect=direct effect + sum of indirect effects
		forvalues i=1(1)4 {
		matrix khbstats[`i',5]=khbstats[1,4]+khbstats[5,3]
		}

		*the contribution percentage of each mediator to the indirect eﬀect, sums up to 100
		forvalues i=1(1)4 {
		matrix khbstats[`i',6]=khbstats[`i',3]/khbstats[5,3]*100
		}
		matrix khbstats[5,6]=khbstats[1,6]+khbstats[2,6]+khbstats[3,6]+khbstats[4,6]

		*the mediation percentage (contribution percentage of each mediator to the total eﬀect, sums up to the overall confounding percentage=(sum of indirect effects)/(total effect))
		matrix khbstats[5,7]=khbstats[5,3]/khbstats[1,5]*100 //overall confounding percentage
		forvalues i=1(1)4 {
		matrix khbstats[`i',7]=khbstats[`i',6]*khbstats[5,7]/100
		}

		matrix drop erase

		*record results to return 
		*order: (1-4:indir effects, 7-9: mediation percentage) 1.absdaggression, 2.absdgetway, 3.absdfriend, 4.absdintelligence 5.total indirect effects 6.direct effect of absdstrength 7.absdaggression, 8.absdgetway, 9.absdfriend, 10.absdintelligence 11.total
		matrix disentangled =J(1,11,.)
		matrix disentangled = [khbstats[1,3],khbstats[2,3],khbstats[3,3],khbstats[4,3],khbstats[5,3],khbstats[1,4],khbstats[1,7],khbstats[2,7],khbstats[3,7],khbstats[4,7],khbstats[5,7]]
			return scalar dagg = disentangled[1,1]
			return scalar dget = disentangled[1,2]
			return scalar dfri = disentangled[1,3]
			return scalar dint = disentangled[1,4]
			return scalar totalind = disentangled[1,5]
			return scalar direct = disentangled[1,6]
			return scalar pdagg = disentangled[1,7]
			return scalar pdget = disentangled[1,8]
			return scalar pdfri = disentangled[1,9]
			return scalar pdint = disentangled[1,10]
			return scalar ptotalind = disentangled[1,11]
restore
end
}
*Run bootstrap
preserve
simulate dagg=r(dagg) dget=r(dget) dfri=r(dfri) dint=r(dint) totalind=r(totalind) direct=r(direct) pdagg=r(pdagg) pdget=r(pdget) pdfri=r(pdfri) pdint=r(pdint) ptotalind=r(ptotalind), ///
  reps(1000) seed(12345) saving(khbdisentangle, replace): myboot2
restore
		


******************************************************************************************************************************************************************************
******************************************************************************************************************************************************************************

**********************************************************
***********Results for mediation analysis from the replication study (Study 3)***********
**********************************************************

cd "..\"

use data_study3,clear
set more off

cd "./results"
quiet{
*----------------------------------------------------------------*
* Recode ratings *
*----------------------------------------------------------------*

*summary stats for raters
*7=extremely unlikely, 1= extremely likely to be prone to aggression
sum rating*

** recode so that higher values mean that characteristic is more likely to fit person in picture**
gen aggressionR=0 if rating<.

replace aggressionR=14 if rating==1
replace aggressionR=13 if rating==2
replace aggressionR=12 if rating==3
replace aggressionR=11 if rating==4
replace aggressionR=10 if rating==5
replace aggressionR=9 if rating==6
replace aggressionR=8 if rating==7
replace aggressionR=aggressionR-7 if aggressionR<.


zscore aggressionR if aggressionR<.


*----------------------------------------------------------------*
* Prepare data for regressions*
*----------------------------------------------------------------*
*Note: in the following double precision is important so that c== actually works, because RHS varable is double precision, whereas standard generate creates single precision and rounding causes == to fail
set type double

*----------------------------------------------------------------*
**Store personality ratings for weaker and stronger opponent in pair, respectively ** 
*----------------------------------------------------------------*

by gameid, sort: gen stronger_aggression = z_aggressionR if strongest==1  
by gameid, sort: gen weaker_aggression = z_aggressionR if strongest==0
by gameid, sort: egen temp = total(stronger_aggression)
replace stronger_aggression = temp
drop temp
by gameid, sort: egen temp = total(weaker_aggression)
replace weaker_aggression = temp
drop temp


*----------------------------------------------------------------*
**Make difference in personality ratings ** (absolute number)
*----------------------------------------------------------------*

by gameid, sort: egen c = min(z_aggressionR) 
by gameid, sort: egen d = max(z_aggressionR)
gen mostaggressive = 0 if z_aggressionR<.
replace mostaggressive = 1 if z_aggressionR==d
by gameid, sort: gen winneraggressive = z_aggressionR if winner==1 & draw!=1
gen mostaggressivewins = 0 if draw!=1
replace mostaggressivewins = 1 if (z_aggressionR==d & winner==1 & draw!=1 )
by gameid, sort: egen temp = max(mostaggressivewins) if draw!=1
replace mostaggressivewins = temp if winner==0 & draw!=1
drop temp
gen absdaggression= abs(d-c)
drop c
drop d


*differences in mediator variables between stronger and weaker contestant*
gen swdiffaggression=stronger_aggression-weaker_aggression

*generate a unique identifier =1 for exactly one line per game (ie per gameid), allows to include each game once only in analyses
by gameid, sort: gen gameno = 1 if _n==1

}


***************************************************************
***********Figure 5: (complete) Mediation***********
***************************************************************
*findit combomarginsplot
which combomarginsplot
preserve
keep if contestantincluded==1
vce2way regress absdaggression absdstrength per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 per11 per12 per13, cluster(id1 id2) 
predict double stragg_res, residuals

quietly: vce2way nbreg finalcosts absdstrength absdaggression per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 per11 per12 per13, cluster(id1 id2)
quietly: margins, at(absdstrength=(0(0.1)3.9)) atmeans saving(fileM1, replace)

quietly: vce2way nbreg finalcosts absdstrength stragg_res per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 per11 per12 per13, cluster(id1 id2)
quietly: margins, at(absdstrength=(0(0.1)3.9)) atmeans saving(fileM2, replace)
combomarginsplot fileM2 fileM1 , labels("Total effect" "Direct effect")  file1opts(recast(line) lw(vthick) lcolor(gray) lpattern(dash)) file2opts(recast(line) lw(vthick) lcolor(black)) noci legend(colfirst) legend(size(medium) ring(0) position(2)) ytitle("Predicted duration in seconds") xtitle("Absolute difference in strength") title("") xscale(titlegap(3))  scheme(lean1) yscale(range(0 60) titlegap(2)) ylabel(0(20)60) saving(filefig4_4, replace)

restore



*Combine figures from studies 3 & 4
gr combine filefig4_3.gph filefig4_4.gph, ///
	col(2) iscale(.7) xcommon scheme(lean1) altshrink  title("       Study 4                                                                        Replication Study(Study 3)", size(small)) 
graph export figure5.pdf, fontface(Arial) replace


********************************************************************************************************************
***** Table S10: (Study 4) Correlations of average personality feature ratings among each other and with strength (standardized) *****
********************************************************************************************************************
pwcorr  z_aggressionR  Formidability_z if gameno==1 , sig obs





********************************************************************************************************************
*Table S12 Mediation in the replication study (Study 3).*
********************************************************************************************************************

which matsave
* if not installed, use findit matsave and load matsave from http://fmwww.bc.edu/RePEc/bocode/m
** WARNING: If the path "dir\matrices\" does not exist, matsave will result in the error **

*Note: this report results based on the pre-saved bootstrap samples
*To run the bootstrap use the code below (search for "RE-RUN Bootstrap")*

cd "..\"


*Prepare a matrix with coefficients, bootstrap standard errors and bias corrected confidence intervals (levels 90,95,99)*
matrix mediation = J(48,6,.)
*								1	2		3		4		5				6						
matrix colnames mediation =  total direct indirect perc indirectisentangle mediationperc  
*								1		2	3	4	5	 6		7	8		
matrix rownames mediation = aggression se ll90 ul90 ll95 ul95 ll99 ul99 ///
getway se ll90 ul90 ll95 ul95 ll99 ul99 friend se ll90 ul90 ll95 ul95 ll99 ul99 ///
intelligent  se ll90 ul90 ll95 ul95 ll99 ul99 ///
totalindir se ll90 ul90 ll95 ul95 ll99 ul99 direct se ll90 ul90 ll95 ul95 ll99 ul99


*** absdaggression*****
quiet{
	*The estimate in our sample
	preserve
	keep if contestantincluded==1
	
	quietly: vce2way regress absdaggression absdstrength per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 per11 per12 per13, cluster(id1 id2)
	quietly: predict double str_res, residuals

	vce2way nbreg finalcosts absdstrength absdaggression per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 per11 per12 per13, cluster(id1 id2)
	matrix erase = e(b)
	matrix khbstats = J(1,3,.)
	matrix khbstats[1,2]=erase[1,1]

	vce2way nbreg finalcosts absdstrength str_res per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 per11 per12 per13, cluster(id1 id2)
	matrix erase = e(b)
	matrix khbstats[1,1]=erase[1,1]
	matrix khbstats[1,3]=khbstats[1,1]-khbstats[1,2]
	matrix drop erase
	matrix list khbstats
	matrix mediation[1,1]=khbstats
	* WARNING: If the path "dir\matrices\" does not exist, matsave will result in the error
	set more off
	matsave mediation, replace p("matrices") dropall
	restore

	set more off
	matload mediation,  p("matrices")  saving over
	
	*Bootstrap se and biase corrected CIs (recognizes what to compute from "estimate in our sample" template above)
	bstat using khbdagg, stat(khbstats) n(205)  level(90)
	estat bootstrap, all 
	
	matrix mediation[2,1]=e(se)
	matrix mediation[3,1]=e(ci_bc)
		
	bstat using khbdagg, stat(khbstats) n(205)  level(95)
	estat bootstrap, all 
	matrix mediation[5,1]=e(ci_bc)
	
	bstat using khbdagg, stat(khbstats) n(205)  level(99)
	estat bootstrap, all 
	matrix mediation[7,1]=e(ci_bc)
}


*** mediation precentage ***
	matrix mediation[1,4]=mediation[1,3]/mediation[1,1]*100

*note: missing values in matrix are -999		
matrix list mediation


*****************************************************************************************************************************
** RE-RUN Bootstrap **
** Run separately to generate bootstrap sample data for table **
*****************************************************************************************************************************


** Write bootstrap samples into results folder (ie don't overwrite original data)**
cd "./results"

*** absdaggression*****
quiet{
capture program drop myboot2
program define myboot2, rclass
	preserve 
		tsset participantid period
		*sample from participant panel observations to "fill" sessions with resampling
		bsample, cluster(participantid) idcluster(newparticid) strata(session)
		*assign new gameid for duplicate games
		by gameid, sort: gen seqrepetition = _n
		egen newgameid=group(seqrepetition gameid)
		*(Not needed for bootstrap) store in each row both subject ids for each game, for clustering
		by newgameid, sort: egen newid1=min(newparticid)
		by newgameid, sort: egen newid2=max(newparticid)
		*generate dummy for which contests to include: drop 205 rows to obtain 205 contests, balanced across periods
		gen randno=runiform()
		sort period randno
		by period: gen pno=_n
		by period: gen pmax=_N/2
		drop if pno>pmax
		drop pno pmax
 
		 regress absdaggression absdstrength per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 per11 per12 per13
		predict double bootstragg_res, residuals

		nbreg finalcosts absdstrength absdaggression per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 per11 per12 per13
		matrix erase = e(b)
		matrix khbstats = J(1,3,.)
		matrix khbstats[1,2]=erase[1,1]

		nbreg finalcosts absdstrength bootstragg_res per1 per2 per3 per4 per5 per6 per7 per8 per9 per10 per11 per12 per13
		matrix erase = e(b)
		matrix khbstats[1,1]=erase[1,1]
		matrix khbstats[1,3]=khbstats[1,1]-khbstats[1,2]
		matrix drop erase
		matrix erase = e(b)
		return scalar khbred = khbstats[1,1]
		return scalar khbfull = khbstats[1,2]
		return scalar khbdif = khbstats[1,3]
	restore
end
}
*Run bootstrap
preserve
	simulate khbred=r(khbred) khbfull=r(khbfull) khbdif=r(khbdif), ///
	reps(1000) seed(12345) saving(khbdagg, replace): myboot2
restore



